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Abstract 

This manual describes a set of utilities developed for Lattice QCD computations. 
They are collectively called qcdutils. They are comprised of a set of Python pro- 
grams each of them with a specific function: download gauge ensembles from the pub- 
lic NERSC repository, convert between formats, split files by time-slices, compile and 
run physics algorithms, generate visualizations in the form of VTK files, convert the 
visualizations into images, perform bootstrap analysis of results, fit the results of the 
analysis, and plot those results. These tools implement the typical workflow of most 
Lattice QCD computations and automate it by enforcing filename conventions: the 
output of one tool is understood by the next tool in the workflow. This manual is 
organized as a series of autonomous recipes which can be combined together. 
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1 Introduction 



In this manual we provide a description of the following tools: 

• qcdutils_get .py: a program to download gauge configurations form the NERSC Gauge 
Connection archive [1] and convert them from one format to another, including to 
ILDG [2] and FcrmiQCD formats [3, 4]. 

• qcdutils_run.py: a program to download, compile and run various parallel Physics 
algorithms (for example compute the average plaquette, the topological charge density, 
two and three points correlation functions). qcdutils_riaii is a proxy for FermiQCD. 
Most of the FermiQCD algorithms and examples generate files that are suitable for 
visuafization (VTK files [5]) 

• qcdutils.vis .py: a program to manipulate the VTK files generated by qcdutils_run 

which can be used to split VTK files into components, interpolate them, and generate 
3D contour plots as JPEG images. This program uses metaprogramming to write a 
Visit [6] script and runs it in background. 

• qcdutils.vtk.py: a program that converts a VTK file into a web page (HTML) which 
displays iso-surfaces computed from the VTK file. The generated files can be visualized 
in any browser and allows interactive rotation of the visualization. This program is a 
based on the "processing.js" hbrary [7] . 

• qcdutils_boot . py: a tools for performing bootstrap analysis of the output of qcdutils_run 
and other QCD Software. It computes autocorrelations, moving averages, and distri- 
butions. 

• qcdutils plot.py: a tool to plot results from qcdutils boot. 

• qcdutils_f it .py: a tool to fit results from qcdutils_boot .py. 

As the .py extension implied, these programs are written in Python [8] (2.7 version recom- 
mended) . 

Together these tools allow automation of the workfiow of most Lattice QCD computations 
from downloading data to computing scientific results, plots, and visualizations. 

Notice that each of the utilities has its own help page which you an access using the -h 
command line option. The output for each is reported in the Appendix. 

The data downloaded by qcdutils_get can be read by qcdutils_rmi which executes the 
physics algorithms implemented in C-|— 1-. The output can be VTK files manipulated by 
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qcdutils_vis and and then transformed into images and movies by Visit, or they can be 
tabulated data that require bootstrap analysis. This is done by qcdutils_boot. The output 
of the latter plotted by qcdutils.plot and can be fitted with qcdutils_f it. 

These files enforce a workflow by following the flle naming conventions described in the 
Appendix but, they do not strictly depend on each other. For example qcdutils.boot can 
be used to analyze the output of any of your own physics simulations even if you do not use 
qcdutils run. 

Here is an overview of the workflow: 

qcdutils-g^ei run > boot ^plot 




vtk 



This manual is not designed to be complete or exhaustive because our tools are in continuous 

development and new features are added every day. Yet is designed to provide enough 
examples to allow you to explore further. Our analysis and visualizations are created on 
sample data and aimed exclusively at explaining how to use the tools. 

Our hope is that these tools will be useful to practitioners in the field and specifically to 
graduate students new to the field of Lattice QCD and looking to jumpstart their research 
projects. 

These tools can also be used to automate the workfiow of analyzing gauge configurations in 
real time in order to obtain and display preliminary results. 

Some of the tools described here find more general application than Lattice QCD and can 
be utihzed in other scientific areas. 

1.1 Resources 

qcdutils can be downloaded from: 
http : //code . google . com/p/qcdutils 
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More information FermiQCD code used by qcdutils_rian is available from refs. [3, 4] and the 

web page: 

http : / /f ermiqcd . net 

More examples of visualizations and links do additional code and examples can be found at: 
http : / /latt iceqcd . org 



1.2 Getting the tools 

There are two ways to get the tools described in here. The easiest way to get qcdutils is to 
use Mercurial: 

Install mercurial from: 



and download qcdutils from the googlecode repository 

h 1 1 p : / / c d e. ^MMMMlM^cmsL/ p/acdutils/source/>--""°°/ ^ 

using the following commands: 

hg clone https://qcdutils.googlecode.com/hg/ qcdutils 



The command creates a folder called "qcdutils" and download the latest source files in there. 



You can also download individual files using wget (default on Linux systems) or curl (default 
on mac systems): 



wget 


http 


/ / qcdutils 


googlecode 


com/hg/ qcdutils 


-get . 


py 


wget 


http 


//qcdutils 


googlecode 


c om / hg / q c du t i 1 s 


_run . 


py 


wget 


http 


//qcdutils 


googlecode 


com/hg/ qcdutils 


_vi s . 


py 


wget 


http 


// qcdutils 


googlecode 


com/hg/ qcdutils 


_vtk . 


py 


wget 


http 


//qcdutils 


googlecode 


com/hg/ qcdutils 


_boot 


•py 


wget 


http 


//qcdutils 


googlecode 


com/hg/ qcdutils 


_plot 


•py 


wget 


http 


//qcdutils 


googlecode 


com/hg/ qcdutils 


_f it . 


py 
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1.3 Dependencies 



These files do not depend on each other so you can download only those that you need. 
qcdutils_nan is special because it is a Python interface to the FermiQCD library. As it is 
explained later, when executed, it downloads and compiles FermiQCD. It assumes you have 
g++ installed. 

qcdutils_f it.py and qcdutils_plot .py requires the Python numpy and matplotlib installed. 
All the file require Python 2.x (possibly 2.7) and do not work with Python 3.x. 

1.4 License 

qcdutils are released under the GPLv2 license. 

1.5 Acknowledgments 
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and code available to the public, and for a long-lasting collaboration. We thank David Skin- 
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instantons. We thank Chris Maynard for useful discussions about ILDG. We thanks Simon 
Catterall, Yannick Meurice, Jonathan Flynn, and all those that over time have submitted 
patches for FermiQCD thus contributing to make it better. We also thank all of those who 
have used and who still use FermiQCD, thus providing the motivation for continuing this 
work. We thank the graduate students that over time have helped with coding, testing, and 
documentation: Yaoqian Zhong, Brian Schinazi, Nate Wilson, Vincent Harvey, and Chris 
Baron. 

This work was funded by Department of Energy grant DEFC02-06ER41441 and by National 
Science Foundation grant 0970137. 

2 Accessing public data with qcdutils _get . py 
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2.1 Searching data on the NERSC Gauge Connection 



The Gauge Connection [1] is a repository of Lattice QCD data, primarily but not limited to 
gauge ensembles, hosted by the National Energy Research Science Center (NERSC) on their 
High Performance Storage System (HPSS). At the time of writing the Gauge Connection 
hosts 16 Terabytes of data and makes it publicly accessible to researchers worldwide. 

The new Gauge Connection site consists of a set of dynamic web pages in hierarchical struc- 
tures that closely mimics the folder structure in the HPSS FTP server. Each folder corre- 
sponds to a web page. The web page provide a description of the folder content, in the form 
of an editable wiki, comments about the content, links to sub-folder and links to files con- 
tained in the folder. Since folders may contain thousands of files, files with similar filenames 
are grouped together into filename patterns. For example all files with the same name but 
different extension or similar names differing only for a numerical value are grouped together. 
Pages are tagged and can be searched by tag. Users can search for files by browsing the folder 
structure, searching by tags and can download individual files or all files matching a pattern. 

You do not need an account to login and you can use your OpenID account, for example a 
Google email account. You do not need to login to search but you need to login to download. 
From now on we assume you are logged in into the Gauge Connection. 

Fig 1 (left) is a screenshot of the main Gauge Connection site. Each gauge ensemble is stored 
in a folder which is represented by a dynamic web page and tagged. You can search these 
pages by tag, as shown in Fig. 1 (right). 



povirered by web2py - @201 1 



C7 rt O qcd.nersc.gov/n^rsc/defajlt^search 




Ths Gauge Connection @ NER5C.gov 





Pages 


o 


Search by Tags 


beta 


o 


flavor 


o 


mass 


o 


source 


o 


volume 


o 





power-ed by web2py - @2011 



Figure 1: Main NERSC Gauge Connection web site (left) and search by tag feature (right). 
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Fig. 2 (left) shows statistical information about tags. 

Each tag has the form "type/value" where the tag type can be: 

• source: the name of the organization who donated the data, value can be, for example 
MILC [9]. 

• flavor, the flavor content of the data, value can be "0" for quenched data, "2" for two 
flavor unquenched data, "2+1" for three flavor unquenched where two quarks have one 
mass and 1 quark has another mass. 

• kappa: the k, value 

• mass: the quark mass 

We have also processed many of the ensambes using some of the tools described here and 
generated animations of the topological charge densities. This is shown in fig. 2 (right). 




The Gauge Connection @ NERSC.go' 



Search by Volume 

V4lumfl/1 2'*3x12 



volunW16'*3xie 



volume^ 6^3x48 



volume/20'^3x2'0 



volume'20'*3x4a 



volume/20^3x64 



Th« Gauge Connection @ NER5C.gov 



CU/D001/u_CU_0001.S400 



CU/0004Ai_C U_0004.SOOO 



IVliLC/1241/u_IWILC_1241.2250 



l\JliLC/1242/u_l\illLC_1242.2450 



IVliLC/1243/u_l\fllLC_1243.2400 



MlLC/1 244/u_MILC_1 244.2450 



powered by web2py - @201 1 



powered by web2pv- @2011 



Figure 2: A page showing statitical information (left) and list of visualzations (right). 
A screenshot of a folder page is shown in fig. 3 (left) 

You can see a description, a list of tags, list of file patterns in the folder, and comments. The 
comments are only visible to logged in users. The login link is at top left of the page. 

A screenshot of a page listing the files in an ensemble is shown in fig. 3 (right) 

You can download an individual files by clicking on the file. 
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X ^^^^^^^^^^^^^^^^^^^^^^^^^^^ 






erst.goWn«rsc/d«rajlt/raot/confTss/MILC/1648f21be503ni049Zni[)S2 






The Gauge Connection @ NER5C.gov 





conf igs/MILC/1 648f21 b6503m0492m082/ 

The requested URL/canfigs/MILC/1 S48f21 b6503m0492mOB2/was not found on this server. 
Apache/2.2.3 (CeritOS) Server at qcd.nersc.gov Port 80 



^ volume/1 6"3x48 flavor/2+1 beta/e.503 mass/0.00492 

I mass/0.0€S2 ^ourcflymilc 



^J.-^ a a O qcd.nersc.goWn«rsc/d«rajlt/files/confl9s/f^lLC/lS48f21bGS[>3niO«2mOS2/u_MILC_lie4Bf?lb„. | 



The Gauge Connection @ NERSC.g 



conf igs/MILC/1 648f21 b6503m0492m082/ 

(400 files, 15.0 GB) 

[ folder I I download tool download link 

(Warning: the files are on tapg and it may take some time to retrieve a file when you dick on a link below. We 
suggest you use the script above instead.) 





u_MILC_l164£f21 b6503m0492m0e2.nnnnn (2010, 400 files. 15.0 GB) 




o 




Your comment: 








B 

u 


_MILC_l1fi4£f21b6S03m0492m0S2.1000(2010-Apr-22, 37.0 MB) 


o 




_IWiLC_l1648f21b6503m0492mOB2.1004 (2010-Apr-22, 37.0 MB) 


o 




.MILC_l1648f21be503m0492m082.1008 (20ia-Apr-22, 37.0 MB) 


o 


u 


_MILC_l1G48f21be503m0492m082.1012 (20ia-Apr-22,37.0MB) 


o 





powered by web2py - @201 1 



powered by web2py - @201 1 



Figure 3: Folder page (left) and page listing all files in an ensamble (right). 



To download files in batch you need to first download qcdutils_get .py. The web page 
page above includes a link [download too/] that explains where to get and how to use 
qcdutils_get .py. We suggest you first read the rest of this section but also read the linked 
instructions which may be more updated. The page also contains a link [download link] which 
is used to reference the data for later download. 

2.2 Downloading data from NERSC 

Here we assume you want to download the 400 MILC gauge configurations of size 16^ x 48 
and (3 = 6.503 computed using 2+1 quarks of mass respectively 0.00492 (light) and 0.082 
(heavy). These files can be found at 

1 http : //qcd . nersc . gov/ners c / def ault /root / conf igs /MILC / 1 648 f 2 Ib6503m0492m082 

where you should notice the folder name 

1 1648f 21b6503m0492m082 

It follows the MILC filename convention 

1 [time] [space] b [beta] m [mass] m [mass] 

and the values of [beta] and [mass] omit the decimal point. 
The above page links the pattern page: 
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1 http://qcd.nersc.gov/nersc/default/files/configs/MILC/1648 
f 21b65 03m0492ni082/u_MlLC_11648f 21b6503m0492m082 . nnnnn 



where 



is the filename pattern and nnnnn is just a wildcard for the gauge configuration numbers in 
the ensemble. 

The page contains a [download link] to a document in JSON format listing all files in the folder 
matching the pattern and additional meta-data about each file. You do not need to open 
this document. All you need to do is copy the link address and pass it to qcdutils_get .py 
as a command line argument. The program opens the URL, download the list, loop over the 
files in the ensemble, and download them one by one. 

Copy the download link to clipboard and it looks something like this: 

1 http://qcd.nersc.gov/nersc/api/files.json/. . ./ 

u MILC 11648f21b6503m0492m082 .nnnnn . 

We have shortened the full actual path using This URL is a personal token and different 

users get different URLs for the same data. This allows the server to monitor usage and 
expire an URL in case of indiscriminate downloads from one user without affecting other 
users. 

To download all data referenced by this link you simply paste the download link after a call 
to qcdutils_get: 

1 python qcdutils_get . py [download link2__^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 

qcdutils_get .py performs the following operations: 

Before downlaod, qcdutils creates a folder with the same name as the ensemble: 
1 u_MILC_11648f 21b6503m0492m082 .nnnnn/ 

and then download all files in there. The files retain the original file name. 

qcdutils also creates a file called "qcdutils. catalog" where it keeps track of successful down- 
loads. This allows automatic resume on restart: if your download is interrupted, for any 
reason (for example network problem or server crash), you can re-issue the download com- 
mand and it resumes where it stopped, qcdutils does not download again files that were 
already downloaded and are currently present on your system. 
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qcdutils can check if a file is complete by checking its size. Data integrity during transmission 
is guaranteed by the TCP protocol. It is still possible that data is corrupted at the source 
or locally after download (for example due to a bad disk sector). If a file is found to be 
corrupted simply delete it, run qcdutils_get again, and it downloads it again. 

Notice that most of the files stored and served by the Gauge Connection are in either the 
NERSC 3x3 or the NERSC 3x2 file format, described later. If your program can read them, 
you do not need any conversion. Yet it is likely you need to convert them and this is the 
subject of the rest of the section. 



2.3 Testing download 

If you encounter any problem downloading real data you can try download a single small 

demo gauge configuration: 

1 python qcd utils_get . p y h ttp : //qcd . nersc . org /n ersc/api/ files/ demo 



It creates a folder called demo and download a single file 



2.4 Converting to ILDG format (.ildg) 

The qcdutils_get.py can also auto-detect and convert file formats. It can input NERSC3x2, 
NERSC3x3, MILC, UKQCD, ILDG, SciDAC, FermiQCD and it can output ILDG and Fer- 
miQCD formats. Other output formats may be supported in the future if this becomes 
necessary, qcdutils converts files using the following syntax 

Here [source] can be a download link, a glob pattern such as "demo/*", or an individual file, 
[target-format] is one of the following: 

• ildg converts a gauge configuration to ILDG 

• mdp converts a gauge configuration to the FermiQCD format 

• slice. mdp converts a gauge configuration (for example 12^ x 48 ) into multiple con- 
figuration files, one for each time slice (for example 1 x 12^ ), in the FermiQCD file 
format. 
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• prop like mdp but converts file propagators from Scidac-ILDG format into the Fer- 
miQCD format. 

• slice.prop.mdp like prop. mdp but converts a Scidac-ILDG propagator into FermiQCD 
time-slice files. 

Most gauge configuration files are very large and require physics algorithms to run in parallel. 
Yet some algorithms, specifically some visualization ones, can work on individual time-slices. 
slice.mdp and slice.prop.mdp allow you to break large files into time-slices for this purpose. 

Here is an example to convert to ILDG format: 

python qcdutils_get .py http ://qcd.nersc .gov/nersc/api/files/demo 
python qcdutils.get . py -c ildg demo/* 

Or in one single line: 

If the source file is "demo/demo. nersc", the converted file has the ".ildg" postfix appended 
and be called "demo/demo.nersc.ildg" . The original file is not deleted. These are the output 
folder /files: 

demo / 

demo/ demo . nerc 

demo /demo .nersc . ildg 



By default qcdutils preserves the precision of the input data, but can specify the precision 
of the target gauge configuration using the -4 fiag for single precision and the -8 fiag for 
double precision. The input precision is automatically detected. For example: 



2.5 Using the catalog file 

The file "qcdutils. catalog" is only used internally by qcdutils and should not be deleted or 
else it loses track of completed downloads and may perform them again unnecessarily. 

If can pass a qcdutils. catalog to qcdutils_get you get a report about the downloaded files, 
python qcdutils_get . py demo/qcdutils . catalog 
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Notice that output files are never overwritten so make sure you delete the old one if you want 
to create new ones. 



2.6 Converting gauge configurations to the FermiQCD format (.mdp) 

This option works similarly to the previous section: 

python qcdutils_get . py http://qcd.nersc.gov/nersc/api/files/demo 



or in one line 

^mibA^ acdrt.il^ crpt. ..^-^.c r.Hn 1^■^■tD://acd.a^x.a^.W^fl^/a1^i/files/dea^ 



It creates 

demo 

demo/ demo . nersc 
demo/demo . nersc . mdp 
demo/qcdutils . catalog . db 

As in the previous case you can specify the precision of the converted file using -4 or -8. 

Notice you cannot specify the endianness. The FermiQCD format (.mdp) uses LITTLE 
endianness by convention because that is the format used internally by x386 compatible 
architectures. 

2.7 Splitting gauge configurations into time-slices 

Often we need to break a single gauge configuration with T time-slices into T gauge config- 
urations with 1 time-slice each. You can do it using the slice. mdp output file format: 

python qcdut il s_get .py http ://qcd. nersc .gov/nersc/api/files/ demo 



The first line creates the following files: 
demo / 

demo /demo .nersc 




i 



HI 



while the second line creates: 
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demo / demo 


. nersc 


. tOOOl , 


. mdp 


demo / demo 


. nersc 


. t0002 , 


. mdp 


demo/demo 


. nersc 


. t0003 , 


. mdp 











Four files because "demo. nersc" contains 4 timeslices. 



2.8 Splitting ILDG progators into timeslices 

We can play the same trick with propagators. While for gauge configurations qcdutils can 
read multiple file formats, for input propagators qcd can only read FermiQCD and SciDAC 
propagators. 

Given a file "propagator.scidac" we can convert it into FermiQCD format: 



which creates 

propagator . scidac . prop .m^ 



or split into time-slices 




This creates 



propagator . 


, scidac , 


. toooo 


. prop 


. mdp 


propagator . 


. scidac , 


. tOOOl 


■ prop 


. mdp 


propagator . 


, scidac , 


. t0002 


. prop 


. mdp 



In this case you can specify the target precision. 



3 Details about file formats 



In this section we show simplified code snippets that should help you understand the different 
file formats used in Lattice QCD. They are very similar to the actual code implemented in 
qcdutils but simplified for readability. 
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3.1 NERSC file format (3x3) 



To better illustrate each data format we present a minimalist program to store data in the 
corresponding format. 

We assume the input is available through an instance of the following class called data: 

1 class GenericGauge ( obj ect ) : 

2 def u (x , y , z , t , mu) : 

3 # u_ij below are complex numbers 

4 return [ [u_00 , u_0 1 , u_02 ] , 
6 [u_10 , u_ll , u_12] , 

6 iiM— Mil \r n "1 IT ""11 ————^^——11 

In this section (and only in this section) we follow the convention that pL — is X, 1 is Y, 
2 is Z and 3 is T. Everywhere else, in particular in the input parameters of qcdutils_run 
fjL — is T, 1 is X, 2 is y and 3 is Z, which is the FermiQCD convention. 

The following code shows how to read data and write it in the NERSC3x3 format: 

1 NERSC_3x3_HEADER = " " " BEGIN_HEADER 

2 HDR_ VERSION =1.0 

3 DATATYPE = 4D_SU3_GAUGE_3x3 

4 DIMENSI0N_1 = %(NX)i 

5 DIMENSI0N_2 = %(NY)i 

6 DIMENSI0N_3 = X(NZ)i 

7 DIMENSI0N_4 = %(NT)i 

8 CHECKSUM = %( checksum) s 

9 LINK_TRACE = %( linktrace) f 

10 PLAQUETTE = %(plaquette) f 

11 CREATOR = %(creator)s 

12 ARCHIVE_DATE = %( archive_date) s 

13 ENSEMBLE_LABEL = %( label) s 

14 FLOATING_POINT = X(precision) s 

15 ENSEMBLE_ID = % ( ensemb I e_id) s 

16 SEQUENCE_NUMBER = % ( s equence_numb er ) i 

17 BETA = %(beta)f 

18 MASS = %(mass)f 

19 END_ HEADER 

20 """ 

21 def save_3x3_nersc ( f ilename , metadata , data) : 

22 f = open ( f ilename , ' wb ') 

23 f . write (NERSC_3x3_HEADER /. metadata) 

24 nt = metadata [ 'ilfr '] 

25 MHHMagMMMMiiHHM 



16 



26 ny = metadata [ 'iVF '] 

27 nz = metadata [ 'iVZ '] 

28 if metaLddLta^L ' PLOATING_POINT ':\== ' IEEE32 ' : 

29 I couple = ' >2f ' 

30 elif metadatal ' PLOATING_POINT ' IEEE64 ' : 

31 couple = ' >2d ' 

32 else : 

33 raise RuntimeError , "Unknown precision" 

34 for t in range (nt) : 

35 for z in range (nz) : 

36 for y in range (ny) : 

37 for X in range (nz) : 

38 for mu in range (0 , 1 , 2 , 3) : 

39 u = data . u (x , y , z , t , mu) 

40 for i in range (3) : 

41 for j in range (3) 

42 c = u [ i ] [ j ] 

43 re , im = r eal ( c ) , imag ( c ) 

44 i f .^w-rAte ( struct . pack ( counle ■ re . im) 1 



The variable couple determines how to pack in binary the 2 variables re , im using big endi- 
anness (">") in single ("f") or double precision ("d"). For more info read the documentation 
on the Python "struct" package. 

All common file formats used by the community to store QCD gauge configuration require 
two loops: one loop over the lattice sites and one loop over the link directions at each lattice 
site. 

In the NERSC, ILDG and MILC case, the first loop is: 

1 f o r t ... 

2 f or z ... 

3 for y ... 

and the second loop is: 

1 for mu in (X,Y,Z,T) # (0,1,2,3) 



3.2 NERSC file format (3x2) 

The NERSC 3x2 format is more common than NERSC 3x3 and here is how to write it: 
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1 NERSC_3x2_HEADER = " " " BEGIN_HEADER 

2 HDR_ VERSION = 1.0 

3 DATATYPE = 4D_SU3_GAUGE 

4 DIMENSI0N_1 = %(NX)i 

5 DIMENSI0N_2 = %(NY)i 

6 DIMENSI0N_3 = %(NZ)i 

7 DIMENSI0N_4 = %(NT)i 

8 CHECKSUM = %( checksum) s 

9 LINK_TRACE = %(linktrace) f 

10 PLAQUETTE = %(plaquette) f 

11 CREATOR = ^Ccreator^s 

12 ARCHIVE_DATE = X ( ar chiv e_ dat e) s 

13 ENSEMBLE_LABEL = %( label) s 

14 FLOATING_POINT = '/.(precision) s 

15 ENSEMBLE_ID = X ( ens emb I e_ i d) s 

16 SEQUENCE_NUMBER = X ( s equence_numb er ) i 



17 BETA = X(beta)f 

IS /f^SS = X(mass)f 

19 END_ HEADER 

20 """ 

21 def save_3x2_ner s c ( f ilename , met adat a , data) : 

22 f = open (filename , 'iwb ') 

23 f . write (MERSC_3x2_HEADER 7. metadata) 

24 nt = metadata [ 'iVT'] 

25 nx = metadata [ 'il^X '] 

26 ny = metadata [ 'iVy '] 

27 nz = metadata [ 'iVZ '] 

28 if metadata [ 'PL0i4rii\^(;_P0Iil^r']== 'J£:££:52 ': 

29 couple = ' >2f ' 

30 elif metadaLtal ' PLOATING_POINT ' IEEE64 ' : 

31 couple = ' >2d ' 

32 else : 

33 raise RuntimeError , "Unknown precision" 

34 for t in range (nt) : 

35 for z in range (nz) : 

36 for y in range (ny) : 

37 for X in range (nz): 

38 for mu in range (0 , 1 , 2 , 3) : 

39 u = data . u (x , y , z , t , mu) 

40 for i in range (3) : 

41 for j in range (2) # here 

42 c = u[i][j] 

43 re , im = real (c) , imag (c) 

44 f . write ( struct . pack ( couple , re , im) ) 



18 



Notice it differs from NERSC 3x3 by only two lines. One line is in the header: 



instead of 

1 DATATYPE = 4D_SU3_GAUGE_3x3 

and in the line marked by a here. 

This file format does not store all the 3x3 matrices but only the first two rows. The third 
row can be reconstructed when reading the file using the condition that the third row is the 
(complex) vector product of the first two. 

qcdutils can reads and rebuilds the missing rows using this code: 



1 


def reunit arize ( s ) : 








2 


(aire , 


alim , a2re , 


a2im , a3re , 


a3im , 




3 


blre , 


blim , b2re , 


b2im , b3re , 


b3im) = s 




4 


clre = 


a2re*b3re - 


a2im*b3im - 


a3re*b2re + 


a3im*b2im 


5 


dim = 


-(a2re*b3im 


+ a2im*b3re 


- a3re*b2im 


- a3im*b2re) 


6 


c2re = 


a3re*blre - 


a3im*blim - 


alre*b3re + 


al im*b3im 


7 


c2iin = 


-(a3re*blim 


+ a3im*blre 


- alre*b3im 


- alim*b3re) 


8 


c3re = 


alre*b2re - 


alim*b2im - 


a2re*blre + 


a2im*blim 


9 


c3im = 


-(alre*b2im 


+ alim*b2re 


- a2re*blim 


- a2im*blre) 


10 


return 


( aire , alim , 


, a2re , a2im 


, aSre , a3im , 




11 
12 




blre , blim , 


, b2re , b2im 


, bSre , b3im , 





3.3 MILC file format 

The MILC file format is the same as the NERSC 3x3 but the header contains different 
information, it uses a binary format, and the endianness is not spedified (although generally 
large endianness is used). The binary data after the header is the same as NERSC3x3. 

Here is an example of code to write a MILC gauge configuration in Python: 

1 def save_milc ( f ilename , metadata , data) : 

2 f = open ( f ilename , 'wb') 

3 milc_header = ' >i^i64.siii ' 

4 milc_magic = 20103 

5 f.write(struct.pack (milc_header , 

6 milc_magic , 

7 metadata [ 'NX '] , 
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8 




metadata [ 


'NY'I , 


9 




metadata [ 


'NZ '] , 


10 




metadata [ 


'#r'] , 


11 




metadata [ 


' ARCHIVE_DATE •'\ [ : 64] , 


12 




metadata [ 


'ORDER '] , 


13 




metadata [ 


'CHECKS UMl '] , 


14 




metadata [ 


•CHECKSUM2 '] ) ) 


15 


nt 


= metadata [ ' nt '] 




16 


nx 


= metadata [' ncc '] 




17 


ny 


= metadata {. 'ny '1 




18 


nz 


= metadata [ ' nz '] 




19 


if 


metadata [ ' PLOATING_POINT <'\== ' IEEE32 ' : 


20 




couple = ' >2f' 




21 


elif metadata [ ' PLOATING_POINT ' IEEE64 ' : 


22 




couple = '>2d' 




23 


else : 




24 




raise RuntimeError , "Unknown precision" 


25 


for t in range(nt): 




26 




for z in range (nz 


) : 


27 




for y in range ( 


ny) : 


28 




for X in range (nz) : 


29 




for mu in r 


ange (0,1,2,3) : 


30 




u = data. 


u (x , y , z , t , mu) 


31 




for i in 


range (3) : 


32 




for j in range (3) 


33 




c = u 


[i] [j] 


34 




re , im 


= r eal ( c ) , imag ( c ) 











When reading a MILC gauge configuration, qcdutils_get checks for the magic number and 
determines the endianness from the first 4bytes of the header, qcdutils also determines the 
precision from the total file size. 

3.4 FermiQCD file format 

The file format used by FermiQCD is called MDP and files have a ".mdp" extensions. They 
are very similar to the MILC format with these differences: 

• the header has a different format and stores slightly different information 

• the endianness is always little-endian. 

• the inner loop over mu has the same order as the outer loop 
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• it is designed to work for an arbitrary number of dimensions (from ID lattices to lOD 
lattices and have an arbitrary site structure) and the FermiQCD code deal with this 
aspect in an automated way that is explored later. 

For regular QCD {SU (3) matrices per link and 4D lattice) a FermiQCD gauge configuration 
can be generated using the following code: 

1 def save_f ermiqcd_4d_ su3 ( f ilename , met adata , data) : 



2 f = open ( f ilename , 'wb') 

3 nt = metadata [' nt '] 

4 nx = metadata [' na; '] 

5 ny = metadata [' 711/ '] 

6 nz = metadata [' nz '] 

7 header_f ormat = ' <60 s60s60sLil0iii ' 

8 maginc_number = 1325884739 

9 ndim , 

10 if ■S!iet&6L&t&\.'PL0ATING_P0INT''\== •IEEE32': 

11 couple = ' >2f ' 

12 metadata [ '5ir£_5IZ£: '] = 4*9*2*4 

13 elif met&A&t&i' PLOATING_POINT • IEEE64 • : 

14 couple = ' >2d ' 

16 metadata [ 'S'ir£'_S'IZ£ '] = 4*9*2*4 

16 else : 

17 raise RuntimeError , "Unknown precision" 

18 f.writeCstruct.pack ( he ader_f ormat , 

19 'File Type: MDP FIELD', 

20 metadata [ 'FIL£'#j4i>^£' '] , 

21 metadata [ • ARCHIVE_DATE [:60] , 

22 magic_number , ndim , nt , nx , ny , nz , , , , , , , 

23 metadata [ ' SITE_SIZE '} , nt *nx*ny *nz ) ) 

24 for t in range (nt): 

26 for z in range (nz) : 

26 for y in range (ny) : 

27 for X in range (nz): 

28 for mu in range (3 , 2 , 1 , 0) : 

29 u = data . u (x , y , z , t , mu) 

30 for i in range (3) : 

31 for j in range (3) 

32 C = u[i][j] 

33 re , im = real (c) , imag (c) 

34 f . write ( struct . pack ( couple , re , im) ) 



Notice the following: 

• The header is binary but uses a string "File Type: MDP FIELD" to identify the file 



21 



format and version. This allows you to identify the file using an ordinary editor, like 
in the NERSC format. 

• It still uses an integer magic number to allow the reader to check the endianness. 

• It requires an ndim variable which is set to 4 because this manual mostly deals with 4D 
fields. 

• It stores in metadata ['SITE_SIZE'] the number of bytes for each lattice site. For single 
precision this is 4 directions times 9 SU(3) matrix elements times 2 (real+complex) 
times 4 (4 bytes for IEEE32, single precision, float) = 288 bytes. It is 576 bytes for 
double precision. 

3.5 LIME file format 

The file formats described so far store metadata in a header which precedes the binary data. 

The LIME data format is different. It is similar to TAR or MIME as scope. It is designed 
to package multiple files into one file. 

A LIME file is divided into segments (sometime called records in the literature although it 
does not strictly conform to the definition of a record because LIME has nothing to do with 
databases). A segment is comprised of five parts: a magic number, a version number, an 
integer storing the size of the binary data, a segment name, and the binary data. 

The magic number identifies the file as a LIME file and the version number identifies the 
LIME version. This information is repeated for each segment. 

Notice that LIME records do not declare the type of the segments and this has to be inferred 
from the name of the segments. One important caveat of LIME is that some segments contain 
binary data, while some contain ASCII strings such as XML. Segments that contain ASCII 
strings are null-terminated and the terminating zero is counted in the size. Binary segments 
are not null-terminated. This is an important detail when reading the data. 

qcdutils contains a class LIME that can be used to open LIME files and read/write segments 
in or out of order. 

A minimalist implementation of the LIME file format is the following: 

1 class Lime ( ob j ect ) : 

2 def __init__ ( self , f ilename , mode , version = 1): 

3 jmmmmasgtimmmmmiutat^ 
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10 

11 

12 
13 
14 
15 
16 
17 
18 
19 
20 
21 
22 
23 
24 
26 
26 
27 
28 
29 
30 
31 
32 
33 
34 
35 
36 



self 
self 
self 
self 
self 

if mode == 'r 
while True : 
header = self 
if not header 
magic , null 
if magic != 



version = version 
filename = filename 
mode = mode 

file = open ( f ilename , mode ) 
records = [] # [ (name , position , size) ] 
mode == ' rb ' : 



or 



file . read (144) 
break 

size, name = struct . unpack ('/ iigJ2Ss ', header) 
1164413355 : 



raise lOError , "not in LIME format" 
name = name [: name . f ind (' \ ') ] 
position = self . f ile . tell ( ) 

self . records . append ((name .position , size)) 
padding = (8 - (size 8)) 1 8 
self . file . seek( s ize + padding , 1 ) 
def read (self , record) : 

(name , position , size ) = self . records [record] 
self . file . seek (position) 
return (name, self .file, size) 
def __iter__ ( self ) : 

for record in range ( len ( self )) : 
yield self . read (record) 
def write ( self , name , data , size = 
position = self . file . tell () 
header = struct. pack('.'ii ql28s 
self . file . write (header) 
self .file.write(data) 

self . file . write (' \0 '*(8 - (size 8)) I 8) 
self . records . append ( (name ,size , position)) 
def close (self ) : 



# in bytes 



None, chunk = MAXBYTES): 
, self . magic ,self .version ,size , name ) 



The actual implementation in qcdutils is more complex because it performs more checks 
and because it can read and write segments even if they do not fit in RAM, which is not the 
case in the example above. 

Here is an example of usage from Python: 

Open a LIME file for writing 
>>> from qcdutils_get import Lime 



Write two records in it 
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i>>> lime . ■write i ' recordl ', ' 01234567 ') 

2 >>> lime.writeC ' record2 ' , 'other binary data ') 

Close it 



Open the file again for reading: 
i>>> lime = Lime (.' test . lime ',' r ') 



Loop over the segments and print, name size, content: 



1 >>> for name , reader , size in lime: 

2 nriTit (n^^^d^^m^M^re^d(size)) 



3.6 ILDG file format 

The ILDG file format uses LIME to package two segments: 

• One segment contains the metadata marked up in XML. 

• One segment contains the binary data, in the same format as in MILC and NERSC 
3x3. 

The XML markup is specified by ILDG for 4D gauge files. Notice that because the first 
segment refers to the second, many programs that read ILDG expect the metadata segment 
to precede the data segment. 

Here is an example of code to write an ILDG file: 

1 def save_ildg (f ilename , metadata , data , If n) : 

2 lime = Lime ( f ilename , ' wb ') 

3 lime . write ( ' ildg -format ' , " " " 

4 <?xml version = "1.0" encoding = "UTF-8"?> 
6 <il dgFormat > 

6 <version >%( VERSION) s </version > 

7 <field>su3gauge</field> 

8 <precision >%( PRECISION) s </preci si on > 

9 <lx>X(NX)s</lx> 



10 <ly>%(NY)s</ly> 

11 <lz>X(NZ)s</lz> 

12 , <lt>X(NT)s</lt> 
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14 Strip () 7, metadata) 



15 nt = metadata [ 'ilfr '] 

16 nx = metadata [ '#Jr '] 

17 ny = metadata [ 'ilfr '] 

18 nz = metadata [ 'iVJf '] 

19 def writer () : 

20 for t in range (nt) : 

21 for z in range (nz): 

22 for y in range (ny) : 

23 for X in range (nz) : 

24 for mu in range (0 , 1 , 2 , 3) : 

25 u = data . u (x , y , z , t , mu) 

26 for i in range (3) : 

27 for j in range (3) 

28 c = u[i][j] 

29 re , im = real ( c ) , imag ( c ) 

30 yield struct . pack ( couple , re , im) 

31 self .lime.write( ' ildg -binary -data ' , writer) 

32 ^^lf Hm^ . write ( 'ildo-dalfl-LfJf ' .If n) 



Notice that this file takes the same arguments as save_3x3_nersc plus an addition one called 
If n. If n stands for lattice file name. 

The If n is intended to be a Unique Resource Identifier (URI) but it not a Universal Resource 
Locator (URL). The prefix Ifn is not a protocol hke http or ftp. 

3.7 SciDAC file format 

The SciDAC format is used primarily for storing propagators. It uses LIME and it packages 
the following segments: 

• scidac-binary-data: the actual binary data 

• scidac-private-f ile-xml 

• scidac-private-record-xml 

We do not describe it here bccuasc this file type is not used by the tools which arc described 
in this manual. Yet we observer that qcdutils_get can convert this files into FermiQCD 
propagators. 
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4 Running physics algorithms with qcdutils_run . py 



qcdutils_run.py is a program for downloading, compiling, and running FcrmiQCD [3, 4]. 
FermiQCD is a library for parallel Lattice QCD algorithms. The library has been improved 
over time and it now includes algorithms for visualization of Lattice QCD data. You can 
learn more about LatticeQCD from refs. [10, 11, 12]. You can learn more about FermiQCD 
from: 

http : / /f ermiqcd . net 

After you download qcdutils, run the following command: 
1 python qcdutils_run . py -download 

This creates a local folder called "fermiqcd" , download the latest FermiQCD source from the 

google code repository: 

1 http :/ / code . googl e . com/p/fermiqcd 

The source include a file "f ermiqcd. cpp" file, which can parse command line arguments and 
run various physics algorithms, some described in this section. qcdutils_run.py compiles 
this source file and stores the compiled one in: 

Notice the .exe extension is used on all supported platforms. 

qcdutils_run requires g++ and you need to install it separately. 

Now you can run physics algorithms with: 
1 python qcdutils_run . py [options] ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 

qcdutils_run.py internally calls fermiqcd.exe and pass its [options] along. 
You can learn more about the FermiQCD options with 

1 python qcdutils_run . py -h 

The output is reported in the appendix but you are encouraged to run it yourself with the 
latest code. 

qcdutils_ruii.py simply passes its command line arguments to "fermiqcd.exe" which parses 
and calls the corresponding algorithms. Some arguments are special (-download, -compile, 
-mpi, -options, -h) because they are handled by qcdutils_run directly. In particular a call 
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to -options introspects the source of "fermiqcd.cpp" and figures out wliicli arguments are 

supported. 

Notice that FcrmiQCD can do more of what qcdutils_run can access. For example it supports 
staggered fermions (including asqtad), staggered mesons, and domain wall fermions. It can 
do visuahzations using those fields too, but that is not discussed here. 

You can fork "fermiqcd.cpp" and force "qcdutils_run" to use your own source code: 
1 python qcdutils_run . py -compile -source myownf ermiqcd . cpp 



4.1 Running in parallel 

There are two ways to run FermiQCD in parallel with qcdutils_ruii.py. On an SMP machine 
you can simply run with the option -PSIM_NPROCS=<nuinber>. Here is an example that loads 
a gauge configuration and computes the plaquette in parallel using 4 processes: 



python qcdutils_run . py -PSIM_NPR0CS=4 \ 











When running in parallel with -PSIM_NPROCS, FermiQCD uses fork to create the parallel pro- 
cesses and uses named pipes for the message passing. Most PCs and workstations do not allow 
dynamic memory allocation of more then 2GB of contiguous space and this creates problems 
when processing large lattices, even if there is enough total RAM available. -PSIM_NPROCS is 
designed to overcome this limitation. 

FermiQCD with -PSIM_PROCS enables you to run parallel processes on one machine even if 

there is only enough RAM to run one of them at time but not all of them concurrently. This 
is because only one of the parallel processes needs to be loaded in RAM at once and the OS 
can automatically switch between processes by swapping to disk. Communications between 
the parallel processes are also buffered to disk and therefore they work as expected. For 
example: 

• qcdutils_run.py -PSIM_NPR0CS=2 forks two processes (0 and 1) 

• pi is put to sleep and pO is executed 

• If pO sends data to pi the data is stored in a named pipe 

• When pO is completed or attempts to receive data it is put to sleep 

• When pO is put to sleep, pi is loaded in RAM and continues execution. 
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• pi can receive the data sent from pO by reading form the named pipe. 

While this is not very efficient, it does allow to run most algorithms even when there is not 
enought RAM available. The communication patters are implemented in ways that avoid 

deadlocks. 

A better option is to use MPI and this is the preferred option for production runs. If you 
want to use MPI, it must be pre-installed on your system separately. On Debian/Ubuntu 
Linux machines this is done with: 

sudo apt -get install mpich2 
cd ~ 

touch . mpd . conf 
chmod 600 .mpd. conf 



In order to use it from qcdutils_run you need to recompile FermiQCD with MPI: 
python qcdutils_run . py -compile -mpi ^^^^^^^ ^ ■imiim iiium ■■i 



This makes an mpi-based executable for FermiQCD: 



You can run it with 

python qcdutils_run . py -mpi=4 \ 

-gauge : start=load : 1 o ad = demo /demo . nersc . mdp -plaquette 

Internally it calls mpirun. 



4.2 General syntax 

The main options of "qcdutilsjrun.py" are: 

• -gauge: creates, loads, and saves gauge configurations 

• -plaquette: computes the average plaquette 

• -plaquette_vtk: generates images of the plaquette density 

• -polyakov: computes Polyakov lines 

• -polyakov_vtk: computes images from polyakov lines 
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• -topcharge: computes the total topological charge 

• -topcharge_vtk: generates images of the topological charge density 

• -cool: cools the gauge configurations 

• -cool_vtk: cools the gauge configurations and save images of the topological charge at 
every step 

• -quark: computes a quark propagator (different sources are possible) 

• -pion: computes a pion propagator (and optionally saves images of the pion propagator) 

• -meson: computes a meson propagator (and optionally saves images of the meson prop- 
agator) 

• -current .static: computes a three points correlation function by inserting a light-light 
between two heavy-light meson operators (and optionally saves images of the current 
density) 

• -4quark: computes all possible contractions of a 4-quark operator between two light 
mesons. 

Each option takes optional attributes in the form :naine=value. All attributes have default 
values. The -pion, -meson and current_static operators take an optional :vtk=true argu- 
ment needed to save the VTK files for visualization. 

Multplc options can be listed and executed together in one run. Although we recommend 
separating the following operations in different runs: 

• Generate gauge configurations, 

• Compute propagators on each gauge configuration. 

• Measure opeartors by reading previously computed gauge configurations and propaga- 
tors. 

The code described here should be considered and example and other cases can be dealt with 
by modifying the provided examples. 

4.3 Creating a cold or hot gauge configuration 

You can create a cold gauge configuration with the following command 
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The -gauge option sets the gauge parameters of FermiQCD. The option is followed by pa- 
rameters separated by a colon. All parameters have default values. 

qcdutils_run.py creates a cold gauge configuration with volume nt=16:nx=4:ny=4:nz=4, 
SU (Nc) with nc=3, and saves it with the name "cold.mdp" . 

The order of the parameters is not important. All parameters have default values. The 
output lists all parameters which are used. 

You can also run 

1 python qcdutils^run . py -gauge : start =hot : nt =16 : nx=4 

to generate a "hot.mdp" gauge configuration. Notice nc=3 is the default. 

4.4 Loading a gauge configuration 

The start attribute of the -gauge option takes four possible values: 

• cold: makes a cold gauge configuration 

• hot: makes a hot gauge configuration 

• instantons: makes a cold configuration containing one instanton and an, optionally, 
one anti-instanton at given positions. 

• load: loads one or more gauge configurations (if more then one, it loops over them) 

When not set, start defaults to load, and FermiQCD expects to load input gauge configu- 
rations. 

In this case, the load attribute of the -gauge option specifies the pattern of the filenames to 
read. 

You can specify one single gauge configuration by filename or multiple configurations using 

a glob pattern (for example "*.mdp"). 

Here is an example that loads all gauge configurations in the "demo" folder and computes 
their average plaqucttc (-plaquette): 

1 python qcdutils_run . py -gauge : start=load : load=demo/* . mdp -plaquette 
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Similarly if you want to download a stream of NERSC gauge configurations and compute tlie 

average plaquette on each of them you can do: 

python qcdutils_get . py -c mdp -4 http://qcd.nersc.org/nersc/api/files/demo 
python qcdutils_run . py -gauge : load=demo /*. mdp -plaquette > run . log 



When loading gauge configurations there is no need to specify the volume since FermiQCD 
reads that information from the input files. 

If you peek into "fermiqcd/fermiqcd.cpp" you can find code like this: 

if ( argument s .haive (" -plaquette ") ) {. 

mdp << "plaquette = " << average_plaquette (U) << endl ; 



Here arguments, have ("-plaquette") checks that the option is present and average .plaquette (U) 
performs the computation for the input gauge configuration U. mdp is the parallel output 
stream and it double as wrapper object for the MPI communicator. 



4.5 Heatbath Monte Carlo 

Whether you start form a cold, hot or loaded gauge configuration you can generate more by 
using the n attribute. In this example: 

Mtitoix. acdutils.-,rua...py--i^g.ajj.Ee.i.s.fcazJLagold : beta,-4^ii.= 10 : therm = 100 • - 

FermiQCD starts from a cold configuration, and using the Wilson gauge action [13] (default) 
generates n=10 gauge configurations. It perform 100 thermalization steps (therm) starting 
from the cold one and then 5 steps separating the one configuration from the next. 

It saves the gauge configuration files with progressive names: 



cold , 


. mdp 




cold 


. 0000 . 


. mdp 


cold 


. 0001 , 


. mdp 


sold 


. 0099 . 





If you want to change "cold" prefix of numbered filename you can specify the prefix attribute 
of the -gauge option. When this attribute is missing, prefix defaults to the name of the 
starting gauge configuration, i.e. "cold" . 
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When you start from hot or cold, FcrmiQCD generates output files in the current working 
directory. If you start from a loaded file, it generates output files (gauge configurations, 
propagators, vtk files) in the same folder as the input files. 

You can use the optional alg attribute to use an improved action or a SSE2 optimized action. 

Here is the relevant code in "fermiqcd.cpp": 

int nconfigs = arguments . get ( "-gau^ie ", "n ", 0) ; 

forCint n=-l; n<iicoiif igs ; n++) -[ 
if(n>=0) { 

int niter = (n==0) ?ntherm : nsteps ; 
if (gauge_action== "lui Ison ") 

WilsonGaugeAction: :heatbath(U, gauge , niter) ; 
else if ( gauge_act ion== "lui I s on_ i/nproi/e d ") 

ii^MMMW WlilMW i iiii Hi— ^^^^^^— 

Use -options to see which algorithms are available. For example you can declare an improved 
gauge action: 

-gauge :action=wilson_improved:beta=. . . :zeta=. . . :u_s=. . . :u_t=. . . 

where Ut, and Us are the parameters of the improved un-isotropic action defined in ref. [14]. 



4.6 Computing a pion propagator 

We define a pion propagator as 

C2[h] = J](7r(0,O)|7r(+i,x)) (1) 

X 

= E E (0|C(0)7U''(0)#(^'^)7i?a'(^'^)|0) (2) 

X ij^apSp 

= EEi'^^'""(^'^)r (3) 

X i,a 

where 

^^^■'"^(t,x)^(0|{g^"(0),g^'^(t,x)}|0) (4) 

is a quark propagator with source at 0. Here a and h label quark flavours, i and j label color 
indexes, a, /3, 5, p label spin indexes. Notice we used the known identity 
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(0| {r{t,^),q^^{0)} |0) = J2^lS*'^"''(t,^hs^ (5) 

pS 

You can compute C2 using the following syntax: 

python qcdut ils_run . py \ 

-gauge : s t art = cold :beta=4:n = 10:steps=5:therni = 100 \ 

iMIMMWWiMiMWMMMIiBil^ ■MiMMWMMiMMBBMI 

qcdutils_run calls "fermiqcd/fermiqcd.exe" which generates 10 gauge configurations and, 
for each, computes a quark propagators with the given values of k and csw using a fast 
implementation of the clover action (another attribute that can be set) and compute the 
pion progator. 

The -quark option loops over the j,/3 indexes and computes the S^^'°'^{t,x). The -pion 
options loops over the i, a indexes and for every t computes the zero momentum Fourier 
transform in x of eq. 2. 

Notice that by default qcdutils saves all the S components. We can avoid it with save=f alse. 
The pion propagator for each gauge configuration can be found in the output log file. 



The output of qcdutils_ruii in this case is looks like the following. 



C2 [0] 


= 14, 


.4746 


C2 [15] 


= 0, 


. 794981 


C2 [0] 


= 14, 


.4746 


C2 [15] 


= 0, 


.794981 



For each t, C[t] takes a different value on each gauge configuration. 
In some of the following example we rely on the output pattern: 



Later we show how to use vtk=true option to save the progagator as function of x and 
visualize it. We also show tools to automate the analysis of logfiles like "run. log". 

If you peek into "fermiqcd.cpp" you find the following code that computes the pion propa- 
gator: 
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1 


f or ( int 


a=0; a<4; a++) 




2 


for(int i=0; i<nc; i++) { 




3 


psi 


= 0; 




4 


if 


(oii_which_process (U. lattice () ,0,0,0,0) = 


=ME) x . set (0,0,0,0) ; 


5 


psi 


(x , a , i) =1 ; 




6 


psi 


. update () ; 




7 


[. . 


.] 




8 


mul 


_invQ (phi , psi , U , quark , abs_precision , rel 


_precision) ; 


9 


[. . 


.] 




10 


if 


( argument s . have (" -p i on ") ) { 




11 




[. . .] 




12 




f orallsitesandcopies (x) { 




13 




for(int b=0; b<4; b++) 




14 




for(int j=0; j<nc; j++) { 




15 




tmp = real ( phi (x , b , j )* con j (phi (x 


,b, j))) ; 


16 




pion [(x(TIME) -tO + NT)°/,NT] += tmp; 




17 




Q (x) += tmp ; 




18 




} 




19 




> 




20 


> 






21 









Notice the field Q which is used in the next section. It is used for 3D visualizations of the 
propagator. 

4.7 Action and inverters 

You can change the action by setting the action attribute of the -quark option to one of the 
following: clover_f ast, clover_slow, clover_sse2. The first of them is the fastest portable 
implementation. The second is a slower but more readable one. The first two support 
arbitrary SU{Nc) gauge groups while the latter is optimized in assembler for Nc — 3. All of 
them support clover, and un-isotropic actions. The attributes are 



kappa 


K, 


kappa.s 


Kg 


kappa.t 




c_sw 


csw 


c_E 


Ce 


c B 


Cb 
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If separate values for Ks,t are not specified, k is used for both, ce is the coefficient that 
multiphes the electric part of the SW term, cb multiplies the magnetic part, csw defaults to 
0. 



The inverter can be specified using the alg attribute of the -quark option and it can be one 
of the following: bicgstab, minres, bicgstabvtk, minresvtk. The meaning of the first two 
is obvious. The second two perform the extra task of saving the field components and the 
residue at every step of the inversion as a VTK file. 

The -quark option also takes an optional source_type attribute which can be point or wall 
and, if point, a source.point attribute to position the source at zero or the center of the 
lattice. It also takes the optional smear _steps and smear_alpha which are used to smear the 
sink. 

The relevant code in "fermiqcd.cpp" is: 

1 forCint a=0; a<4; a++) 

2 for(int i=0; Knc ; i + + ) { 

3 if ( source_type== ' 'point ' ') { 

4 psi = ; 

6 if ( on_whi ch_pr o cess (U . latt i ce ( ) , to , xO , yO , zO ) ==ME ) { 

6 X . set (to , xO , yO , zO) ; 

7 psi (x , a , i) =1 ; 

8 > 

9 } 

10 [ . . . ] 

11 psi . update () ; 

12 [ . . . ] 

13 if ( argument s . get ( ~ ~ - quark '',''load.'',''false| true '')==' 'true '' ) { 

14 phi . load ( quarkf ilename ) ; 

15 } else { 

16 mul_invQ (phi , psi , U , quark ,abs_precision ,rel_precision) ; 

17 phi . save ( quarkf ilename ) ; 

18 } 

19 [ . . . ] 

20 if (use_propagator ) { 

21 f orallsites (x) { 

22 forspincolor(b,j,nc) { 

23 S(x,a,b,i,j) = phi(x,b,j); 

24 } 

25 y 

26 )■ 



27 
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Notice that in FermiQCD inverters are action agnostic. A call to mul_Q(phi,psi,U, . . .) 
computes = Q\U]ip where Q is the selected action for the type of fermion ip (in this 
document we deal only with wilson type fermions but it works with staggered and domain 
wall too). A call to mul_invQ(phi,psi,U, . . .) computes = Q~^\U\'^ using the same Q and 
the selected inverter. There is no code in the inverter which is action specific. 

4.8 Meson propagators 

Given a meson created by qVq |0), a meson propagator can be defined as follows: 

c^ih] = E(r'°""'(o,o)|r^"^(+t,x)> (6) 

X 

= E E (oic(o)cr-gf(o)gf(t,x)q;,"'=g:''(t,x)|o) (7) 

X ij,a0Sp 

= E E {^''''''-f^)spS*''''"'{t, x) (7^r^°"''"")„^ (8) 

X 

The command to compute an arbitrary meson propagator and reuse the previously computed 
propagators (the code assumes different flavours of degenerate quarks, i.e. same mass): 

1 python qcdutils_run . py \ 

2 -gauge : St art = cold : beta=4 : n = 10 : steps =5 : therin = 100 \ 

3 -quark : kappa=0 . 11 : c_sw=0 . 4 : save=f alse \ 

4 -geson : source_gamma = l : sin k_gamm a = l > j;u jL^. 1 ^^^^^ 

The source-gamma and sink_gamma attributes can be specified according to the following table: 
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source_gcLmnia/ ri nk gamma 




T 

i 


1 
i 


r 



7 


u 


7 


1 
i 


7 




7 


Q 
O 


7 


Do 


7 7 


1 K 

io 




AO 


7^7^ 


OD 


7^7^ 


m 
Ui 


7 7 


uz 


7 7 


03 




12 


7^^ 


13 




23 


7^7^ 



The relevant code in "fermiqcd.cpp" is described here: 



1 if ( arguments . have ( "-meson ") ) { 

2 [...] 

3 Gl = GammaS *par se_gaiiiiiia ( argument s . get (" -mes on ", "s ourc e_ 5 awma ) 

4 G2 = par se_gamma ( argument s . get (" -mes on "," s infc_5a/n/na ",...))* GammaS 
6 f orspincolor (a , i , U . nc ) { 

6 f or spincolor (b , j , U . nc ) { 

7 f orallsites (x) { 

8 sl=s2=0; 

9 for(int c=0;c<4;c++) { 

10 si += S (x , a , c , i , j ) *G2 (c , b) ; 

11 s2 += conj (S (x , c ,b , i , j ) ) *G1 (c , a) ; 

12 } 

13 tmp = abs(sl*s2); 

14 meson [(x(TIME) -tG+NT)%NT] += tmp; 

15 Q(x) += tmp; 

16 } 

17 _ > 



As before we use a scalar field Q for data visualization. 

In this and the other examples the two quarks are degenerate but it is possible to change one 
of the quark propagators by simply replacing it in the code for a different one. We leave it 
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to the reader as an exercise. A next version of "fermiqcd.cpp" will have an option -quark2 
for doing this automatically. 

4.9 Current insertion 

We define it as follows (for two light quarks a, b and one static quark h): 

CcurreuM^ = J] x) | ^ar™*g,(0) | x) > 

X 

= E E ^"(-^' ^)ra^?f (-i, x)g;'fr^r'-*gf (i, x)r,,/i^''(i, x) |0) 

X 

= ^tr(P°"''^Y'5n-^>x)7^r'=™*5(t,x)r*"'=i/t(-i,i,x)) (9) 

X 

Here H is the heavy quark propagator according to Heavy Quark Effective Theory [15] (from 
(-i,x) to (i,x)): 

H{-t,t,x) = ^{1 + -f'')Uo{-t,x)Uoi-t + l,x)...Uoit - l,x) (10) 



You can compute it with 



python qcdutils_run . py \ 

-gauge : start=cold:beta=4:n=10: steps=5: 
-quark : kappa=0 .ll:c_sw=0.4:save=false 


:therin = 100 \ 
\ 






The relevant code in "fermiqcd.cpp" is: 


Gl = parse_gainiiia ( arguments . get ( "-cwrrent. 


.static", " source_ 


gamma " , . 


. .))* 



GanunaS ; 

G2 = par se_gamma ( arguments .get( " -cur rent _static" , " sink_gamma ",...)); 

G3 = GammaS *par se_gamma ( arguments . get ( " - current_ static " , " current _ g amma " 



,...)); 

G4 = G2*(l-Gamma [0] ) /2*G1 ; 
f orallsites (x) 

if (x(TIME) >=0) { 

z.set((NT+2*tO-x(TIME))%NT,x(l) ,x(2) ,x(3)) ; 
f orspincolor (a , i , U . nc ) { 
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9 f or spincolor (b , j , U . nc ) { 

10 sl = s2 = 0; 

11 for(int c=0; c<4; C++) { 

12 si += conj (S (z , c , a , j , i) ) *G3 (c ,b) ; 

13 forCint k = 0; k<U.nc; k++) 

14 s2 += S (x ,b , c , j , k) *G4 (c , a) *conj (Sh (x , i , k) ) ; 

15 } 

16 tmp = abs(sl*s2); 

17 current [(x(TIME)-tO + NT)7.NT] += tmp; 

18 Cj(x) += tmp; 

19 > 

20 y 



Here Sh is the product of links from —t to t along the time direction. 



4.10 Four quark operators 

Instead of inserting a current we can insert a 4-quark operator between two meson operators 
(hght-light): 

C3[ti][t2] = J](r™(-ti,xi)|g„r^g„®g,rBgd|P"''(+t2,x2)> (11) 

Xl X2 

= tr(r— T^-^^l-^i, xi)7'r^-s(-ti, xi))tr(r*"V^^(t2, ^2h'rBS{t2, xO) 

or tr(r— 7^,5t(-ii, xi)7^r^,5(i2, ^2)r'''"'^'S^ {h, X2)7'rB-S(-ii, Xl)) 

The or indicates that there are two possible contrations. FermiQCD computes both of them 
and writes them seperately in the output. 

Here Ta ® Tb is the spin/color structure of the 4-quark operator. We are also ignoring the 
contractions that corresponds to disconnected diagrams. 

We can compute ® Tb for 75 (g) 75 in spin and 1 ® 1 in color (5Ix5l) with: 

1 python qcdutils_run . py \ 

2 -gauge : start=cold : beta=4 : n=10 : steps=5 : therm=100 \ 



3 



In this example, source=l indicates that p*""^^^ — F*'"'^ — 7^. 

This generates the following output, repeated for each of the 10 gauge configurations: 
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1 C3 [0] [0] = 9 . 12242 

2 C3 [0] [1] = .485189 

- ^-^ri.KiL.i^j^ q 1-^-^^-^ , 

Notice the program computes the two contractions of the operator and writes one in C3 and 
one in C3x. 

Instead of source=l you can use any of the operators defined for mesons. 

Instead of 51x51 4-quark operator you can use any the following other operators: 51x51, 
01x01, 11x11, 21x21, 31x31, 051x051, 151x151, 251x251, 351x351, OlIxOlI, 021x021, 031x031, 
121x121, 131x131, 231x231, 5Tx5T, OTxOT, ITxlT, 2Tx2T, 3Tx3T, 05Tx05T, 15Txl5T, 25Tx25T, 
35Tx35T, OlTxOlT, 02Tx02T, 03Tx03T, 12Txl2T, 13Txl3T, 23Tx23T. Here the numerical part 
represents the T^T stucture of the 4-quark operator, the J or T represents its color structure. 
TxT stands for Ea^" ® ^i^^ = and A'' is the SU{3) generator. 



Here is the relevant source code in "fermiqcd.cpp": 



1 


mdp_niat rix G = parse _ gamma ( argument s . ge t ( " quavk " , "source", . 


. . ) ; 


2 


f or spincolor ( a , i , U . nc ) { 




3 


forCint c=0; c<4; C++) 




4 


forCint d=0; d<4; d++) 




5 


if (G(c ,d) !=0) 




6 


f orallsites (x) 




7 
8 


forCint k=0; k<U.nc; k++) 

open_prop [a] [b] [i] [j] [ (x (TIME) -tO + MT) 7.NT] += 




9 


S(x,a,c,i,k)*conj (S(x,b,d,j ,k))*G(c,d) ; 




10 


string op4q = arguments . get ( "-^guarfc ", "operator ",...) ; 




11 


if ( arguments . have ( " -4 quark ") ) {, 




12 


forCint a=0; a<4; a++) 




13 


forCint b=0; b<4; b++) 




14 


forCint c=0; c<4; C++) 




15 


forCint d=0; d<4; d++) { 




16 


mdp_complex gl = GlCb.a); 




17 


mdp_complex g2 = G2Cd,c); 




18 


if Cgl ! =0 && g2 ! =0) 




19 


forCint i=0; KU.nc; i++) 




20 


forCint j=0; j<U.nc; j++) 




21 


if C ! rotate ) { 




22 


c3a+=abs Copen_prop [a] [b] [i] [i] [tls]*gl* 




23 


open_prop [c] [d] [j] [j] [t2s]*g2) ; 




24 


c3b+=abs Copen_prop [c] [b] [j] [i] [tls]*gl* 




25 


open_prop [a] [d] [i] [j] [t2s]*g2) ; 




26 


} else 
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27 [ . . . ] 

28 } 

29 y 

30 } 

31 } 



Notice the two contractions are computed separately. The case rotate==true corresponds to 
the TxT color stucture. 



5 Images and movies with qcdutils vis .py and qcdutils_vtk.py 

In this section we describe how to make 3D visualizations using Visit [6] and how to embed 
visualizations into web pages using "processing.js " [7] . 

Visit is a visualization software developed at Lawrence Livermore National Lab based on 
the VTK toolkit. It provides a GUI which can be used to open the VTK files created by 
FermiQCD (or other scientific program) in interactive mode, but it can also be scripted using 
the Python language. 

"processing.js" is a lightweight javascript library that allows drawing on an HTML canvas 
using the processing language or the javascript language. 

qcdutils uses meta-programming to generate Visit scripts (qcdutils.vis) or processing.js 
scripts (qcdutils_vtk). The former is more flexible and is more appropriate for making high 
resolution images. The latter makes it easy to embed 3D visualizations into web pages. 

Using Visit is intuitive but there are certain tasks which can be repetitive. For example if 
you have multiple VTK files containing topological change density (or any other scalar field), 
you have to determine the optimal threshold values for the contour plots. If you have many 
files you may want to interpolate between them for a smoother visualization, qcdutils.vis 
helps with these tasks. In particular it can: 

• Spht VTK files containing multiple time-slices into separate VTK files, one for each 
shce. 

• Interpolate each couple of consecutive VTK files and make new ones in between. This 
is necessary for smoother visualizations. 

• Compute automatic thresholds values for contour plots. 

• Resample the points by interpolating between the. 
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• Generate Visit scripts which converts VTK files to JPEG format (these script can be 

saved, edited, and reused). 

• Pipe the above operations and run them for multiple files. 

Images generate in this way can be assembled into mpeg4 (or quicktime or avi) movies using 
ffmpeg (an open source tool that is distributed with Visit) but there are other and better 
tools available. We strongly recommend "MPEG Streamclip" . It is much faster, robust, and 
much easier to properly confgure than ffmpeg. 

5.1 About VTK file format 

There are many VTK file formats, qcdutils uses the binary VTK file format described below 
to store scalar fields, usually by timeshces. 

A typical file has the following content: 

1 # vtk DataFile Version 2.0 
■2 f ilename . vtl 

3 BINARY 

4 DATASET STRUCTURED_POINTS 
.-, DIMENSIONS 4 4 4 

6 ORIGIN 

7 SPACING 111 
s POINT_DATA 64 

9 SCALARS scalars_tO float 

10 LOOKUP_TABLE default 

11 [binary data] 

12 SCALARS scalars_tl float 

13 LOOKUP_TABLE default 

14 [binary data] 

15 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 

It consists of an ASCII header declaring the 3D dimensions (4 4 4) and the total number of 
points (4 X 4 X 4 = 64 ). This is following by blocks representing the time-shces. Each block 
as its own ASCII header: 

1 SCALARS scalars_tO float 

2 LOOKUP_TABLE default 

followed by binary data (64 floating point numbers). 



42 



scalars_tO, scalars_tl, etc. are the names of the fields as stored by FermiQCD. When 
time-slices are extracted by qcdutils_vis the slices are renamed as slice. 

Given any VTK file, for example demo.vtk we can visualize it using qcdutils_vis .py using 

the following syntax: 

1 python qcdutils_vis . py -r ' scalars_tO ' -p default demo . vtk 

qcdutils_vis .py generates images in JPEG format. 

Similarly we can visualize by creating an interactive 3D web page: 

If the filename is a glob pattern (*.vtk), both tools loop and process and files matching the 
pattern. 

qcdutils_vtk computes the range of values in the scalar field from the maximum to the 
minimum, -u 0.10 indicates we want an isosurface at 10% form the max and -1 0.90 
indicates we want another isosurface at 90% from the max (10% from the min). It is also 
possible to specify the colors of the iso-surfaces. 

qcdutils_vtk generates HTML files with the same as the input VTK files followed by the 
.html postfix. The isosurfaces are computed by the Python program itself but the repre- 
sentation of the isosurfaces is embedded in the html file, together with the "processing.] s" 
library, and with custom JS code. These files are not static images. You can rotate them in 
the browser using the mouse. 

5.2 Plaquette 

As an example, we want to make a movie of the plaquette as function of the time-slice. We 
follow this workflow: 

• Load a gauge configuration. 

• Gompute the plaquette at each lattice site. 

• Save the plaquette as a VTK file. 

• Split the VTK file into one file per time-slice. 

• Interpolate the timeslices to generate more frames. 

• Generate contour plots for each frame and save them as JPEG files. 
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This can be done in two steps. Step one: 

python qcdut ils_run . py \ 

-gauge : load=demo/demo . nersc . mdp \ 



This command uses FermiQCD to load the gauge configurations. For each of them it com- 
putes the trace of the average plaquette at each lattice site, and generates one VTK file 
contain the 4D scalar for the plaquette. This file is saved in a new file with the same prefix 
as the input but ending in ". plaquette. vtk" . 

Step two: 

PinHHipMHMHPPHHHHHII^HHMH^^ 

It reads all files matching the pattern "demo/*. plaquette. vtk" , extracts all time-slices with 
names matching "*" (all time slices), and interpolates each couple of VTK files by adding 
9 more frames (-i 9), then generates a VTK script that reads each VTK file, resamples it, 
and stores contour plots in JPEG files with consecutive file filenames.. 

The generated script has a unique name which looks like this: 
qcdut ils_vis_2faclb86 -5b86 -42ee -8552 . py 

qcdutils_vis writes and runs the script. It saves it for you in case you want to read and 
modify it. 

When it runs, it loops over all the frames, resamples them, computes the contour plots and 

saves each frame into one JPEG image: 

qcdut ils_vis_2faclb86 -5b86 -42ee-8552_0000.00.jpeg 
qcdut ils_vis_2faclb86 -5b86 -42 ee -8552 _0001 .01 . jpeg 

qcdutils_vis_2f aclb86 -5b86-42ee-8552_0003 .00. jpeg 



Here 0000, 0001, 0002, 0003 are the original frames (timeslices) and the. .01, .02, .09 
are the interpolated ones. 

Notice that the -i 9 option is very important to obtain smooth sequences of images to be 
assembled into movies. 

The option 
is equivalent to 
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1 Mlii'iliMlfilif^'^^^^-*^ " : R.,,....a^l.^^^J...i.^JLL^^.,^^ 



Here Annotation, Resample and Contour are Visit functions. Using -p you can set the at- 
tributes for each functions. 

For example, to remove the bounding box you would replace 
1 AnnotationAttributes [] 



with 



■Mi 



To increase the re-sampling points from 100 to 160 you would replace: 
1 ResampleAttributes [] 

with 
1 ^^^H 

To change the color of the 9th contour to Orange, you would replace: 
1 Contour Attributes [] 



with 



The argument of the <f unction> Attributes [. . .] are Visit attributes and they are described 
in the Visit documentation. 

The relevant page of code in "fermiqcd.cpp" that computes the VTK plaquette is here: 

1 void plaquette_vtk ( gauge_f ieldfe U, string filename) ■[ 

2 mdp_f ield <mdp_real > Q (U . lattice ()) ; 

3 mdp_site x (U . lattice ()) ; 

4 f orallsites (x) if(x(0)==0) { 

5 Q(x)=0; 

6 for(int mu=0; mu<4; mu++) 

7 for(int nu=mu+l; nu<4; nu++) 

8 Q(x)+= real (trace (plaquette (U , x , mu , nu) ) ) ; 

9 > 

10 Q . save_vtk (filename , -1) ; 

11 } 

12 [ . . . ] 

13 if ( argument s . have ( "-p I ague 1 1 e_ tfe ") ) { 
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Figure 4: Visualization of a contour plot for the average plaquette (left) and the intersection 
of the contours with the bounding box (right) 

14 plaquette_vtk(U,newfilename+".plagwette.t/tfc") ; 

15 } 

Notice how the plaquette is computed for each x, summed over mu,nu, stored in a scalar field 
Q(x), and then saved in a file. This strategy can be used to visualize any FermiQCD scalar 
field with minor modifications of the source. 

5.3 Topological charge density 

Similarly to the average plaquette we can make images corresponding to the topological 
charge density. 

To generate the topological change density we need to cool the gauge configurations (-cool) 
and then compute the topological charge (-topcharge_vtk): 



python 


qcdutils_run . py \ 

-gauge : load=demo/demo . nersc . mdp \ 
- cool : steps =20 -topcharge_vtk 






python 


qcdutils_vis . py -s -i 9 -p default 


' demo / * . topcharge . vtk ' 



The relavent code in "fermiqcd.cpp" is below: 
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Figure 5: Visualization of the topological charge density. 

1 if (ar guments . have ("- t op char g e_vtk ") ) { 

2 float tc = t opologi cal_ charge_vtk (U , newf ilename +". t op cftar^e . 'utfc ",- 1) ; 

3 mdp << "topcharge = " << tc << endl ; 

4 > 

topological_charge_vtk is defined in "fermiqcd_topological_charge.h" . The -1 arguments 
indicates we want to save all time slices. The actual code to compute the topological charge 
density is: 

1 void topologi cal_charge (mdp_f ield <f loat > &Q , gauge_field &U) { 



2 compute_em_notrace_f ield (U) ; 

3 mdp_site x (U . lattice ()) ; 

4 f orallsitesandcopies (x) { 

5 C!(x)=0; 

6 forCint i=0; i<U.nc; i++) 

7 forCint j=0; j<U.nc; j++) 

8 Q(x)+=real(U.ein(x,0,l,i,j)*U.em(x,2,3,j,i)- 

9 U.ein(x,0,2,i,j)*U.em(x,l,3,j,i) + 

10 U.em(x,0,3,i,j)*U.em(x,l,2,j,i)); 

11 } 

12 Q . updat e ( ) ; 



13 > 

Here U.em is the eletro-magnetic field computed from U. 
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5.4 Cooling 

Sometimes we may want to see the changes in the topological charge density as the con- 
figuration is cooled down. This requires computing the topological charge density at every 
cooling step. This can be done with the -cool.vtk option: 



python 


qcdutils_ruii . py \ 

-gauge : start=load : load=demo/ demo 

-cool_vtk : cooling=10 > run. log 


. ner sc . mdp \ 






acdutils,^.¥ls...Dv -r ' scalars tO ' 


-i 9 -p default ' demo /*, cool ??. vtk ' .. 



The -cool_vtk option creates VTK files ending in "coolOO.vtk" , "coolOl.vtk" "cool49.vtk" . 
To make a smooth movie we do not break files into time-slices (no -s option) but instead we 
extract the same slice for every file (-r 'scalars_tO). Then we interpolate the frames (-i 

9)- 

The above code generates JPEG images showing different stages of cooling of the data. You 
can see some of the images in fig. 6 

The relevant code in "fermiqcd.cpp" is here: 

void cool_vtk (gauge_f ieldfe U, mdp_args& arguments, string filename) { 
if (arguments . get ( "-cool " , "alg " , "ape ") == "ape ") 

for(int k=0; k<arguments . get ( "-coo l.t^ifc ", "n ", 20) ; k++) { 
ApeSmearing : : smear (U, 

argument s . get (" - co o l_vtk " , "alpha ",0.7) , 
arguments .get("-cooi_vtfc", "steps " ,1) , 
argument s .get( " - m n i m + h " ^ " n nli.ng " ,10) ) ; 
topological_charge_vtk(U,filename+". cool "+ to string (k,2) + ".i;tfc",0); 

} 

else 

mdp . err or _mes sage (" CO o I in^ algorithm not supported"); 



The smearing algorithm is in the "topologicaLcharge_vtk" file: 

class ApeSmearing { 

public: static void smear ( gauge_f ield &U , 

mdp_real alpha=0.7, 
int iterations=20 , 
int cooling_steps=10) { 
gauge_field V (U . latt i ce () , U . nc ) ; 
mdp_site x (U . lattice ()) ; 

for (int iter=0; iter < iterations ; iter++) { 
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Figure 6: Visualization of the topological charge density at different cooling stages. 



9 


cout << "smearing step " << iter << "/" 


<< iterations << endl ; 


10 


V=U; 




11 


forCint mu = 0; inu<4; mu++) { 




12 


f orallsites (x) { 




13 


U (x , mu) = ( 1 . 0- alpha) *V (x , mu) ; 




14 


forCint nu=0; nu<U.ndim; nu++) 




15 


if (nu ! =mu) 




16 


U(x,mu )+=(!. 0-alpha)/6* 




17 


(V(x , nu) *V(x+nu , mu) *herinitian 


(V (x+mu , nu) ) + 
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18 








hermitian (V(x-nu,nu))*V(x-nu 


,mu)*V((x-nu)+mu,nu)) ; 


19 




U(x, 


mu) 


=pro j ect_SU (U (x , mu) , cooling_ 


steps ) ; 


20 




> 








21 


> 










22 


U, 


. update 


() ; 






23 


} 










24 } 












25 } ; 













5.5 Polyakov lines 

A Polyakov line is the trace of the product of the gauge links along the time direction, 
therefore it is a 3D complex field. Here we are interested in the real part only (the image 
part is qualitatively the same). 

We can visualize Polyakov lines using the -polyakov_vtk option: 

1 python qcdutils_run . py \ 

2 -gauge : load = deino/demo . nersc . mdp \ 

3 iMHHHiiiiiittiiHiiii^^ 

which we can convert to images with: 



1 python qcdutils_vis . py \ 




The output is show in fig. 7. 

Here is the relevant code in "fermqcd.cpp": 



1 


void polyakov_vtk ( gauge_f ieldfc U, string filename) { 




2 


int L [3] ; 




3 


L [0]=U. lattice () . size (1) ; 




4 


L[l]=U.lattice() .size(2) ; 




5 


L[2]=U.lattice() .size(3) ; 




6 


mdp_lattice space (3, L, 




7 


def ault _part it i oning <1> , 




8 


torus_topology , 




9 


, 1 , false) ; 




10 


mdp_matrix_field V( space ,U.nc ,U.nc) ; 




11 


mdp_f ield <mdp_real > Q(space,2); 




12 


mdp_site x (U . lattice ()) ; 




13 


mdp_site y(space); 




14 
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Figure 7: Visualizaiton of contour plots for the Polyakov lines. Different colors represent 
positive and negative values of the real part of the Polyakov lines. 

15 int k , mu = , nu = l ; 

16 mdp_complex s = 0; 

17 

18 f orallsites (y) V(y)=l; 

19 forCint t=0; t<L[0]; t+ + ) { 

20 f orallsites (y) { 

21 X. set (t ,y (0) ,y (1) ,y(2) ) ; 

22 V(y)=V(y)*U(x,0) ; 

23 } 

24 } 

25 

26 f orallsites (y) { 

27 indp_ complex z = trace ( V (y ) ) ; 

28 Q (y , 0) =real (z) ; 

29 Q (y , 1) =imag (z) ; 

30 } 

31 Q. save_vtk(filename ,-1,0,0, false) ; 

32 > 

33 [ . . . ] 

34 if (arguments . haLve (" -p ly akov_vtk ") ) { 

35 polyakov_vtk(U,newfilename + ". polyakov . vtk") ; 

36 } 

This code is a little different than the previous one. It creates a 3D lattice called space which 
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is a time projection of the 4D space. While x hves on the full lattice, y leaves only on the 
3D space, g is a scalar field with two components (real and imaginary part of the Polyakov 
lines) which lives in 3D space. 

5.6 Quark propagator 

Given any gauge configuration we can visualize quark propagators in two ways. We can use 
the normal inverter and save the proprgator at the end of the inversion for each source/sink 
spin/color component: 

python qcdutils_run . py \ 

-gauge : start = load : load = d.eino/demo . nersc . mdp \ 

-quark : kappa=0 .135: sour ce_point=center:alg=bicgstab: vtk=true I run .log 

(here using the bicgstab, the Stabilized Bi-Conjugatc Gradient). Alternatively we can use 
a modified inverter which saves the components but also VTK visualization for the field 
components and the residue at each step of the inversion. 

python qcdutils_run . py \ 

-gauge : start=load : load=demo/ demo .nersc. mdp \ 
llWIIinimiinr^^ • kanD.a=..Q.. .13,5 : : a.lg = bi.c.gst||| yj,^ ^Mlliillliiiiffl^MMi 

Fig. 8 shows different components of a quark propagator on a hot and a cold configuration. 
From now on we assume the propagator has been computed and we reuse it. 

5.7 Pion propagator 

In a previous section, computed the zero momentum Fourier transform of the pion propagator. 
Now we want to visualize it for every point in space: 

(5(t,x) = (7r„6(0,O)|7r,6(+t,x)) = J] |5^^'"«(i, x) |' (12) 

This can be done using the vtk=true attribute of the -pion option: 

python qcdutils_run . py \ 

-gauge : start=load : load=demo/demo . nersc . mdp \ 
-quark : kappa=0 . 135 : source_point=center : load=true \ 
-pion : vtk=true > run. log 
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Figure 8: Different components of a quark propagator on a cold gauge configuration (left) and 
on a thermalized gauge configuration (right). From top to bottom, they show the magnitude 
of 5°''*^ (t,x) = ^oo™(0,x)), 5°200(0,x), and ^°=^™(0,x). 

Notice the -quark . . . : load=true which reloads the previous propagator. We can now convert 
the pion VTK visualization into images using qcdutils.vis: 

1 python qcdutils_vtk . py -u 0.01 -1 0.00001 ' demo /*. pion . vtk ' 
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python qcdutils_vis . py -s '*' -i 9 -p default ' demo /* . pion . vtk ' 

In this case the -i 9 option is used to interpolate between time-shces in case the images are 
to be assembled into a movie. 

Examples of images are shown in fig. 9 




Figure 9: Contour plot for a pion propagator. 



5.8 Meson propagators 



A meson propagator is defined similarly to a pion propagator but it has a different gamma 
structure: 



g(t,x) = (r™(o,o)|rr'(+t,x)> (13) 

= ^^^^■•'''^(t,x)(r^^"V)<5p^"^''"''(t,x)(7^P™)^;9 (14) 
We can visualized a Meson propagator using the following code: 



1 


python 


qcdutils_run . py \ 






2 




-gauge : start=load: load=demo/demo 


. ner s c 


. mdp \ 


3 




-quark : kappa=0 . 135 : source_point= 


center 


: load = true \ 


4 




-meson : sour ce_gamma=l : s ink_gamma 


=1 : vtk 


=true > run . log 
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and then process the VTK file as in the pion example. In this case r*""'"'^^ — — 
indicates a verctor meson polarized along the X axis. 

5.9 Current insertions 

We can also visualize the mass density and the charge density of a heavy-light meson by 
inserting an operator {qq and q^^q respectively) in bewteen meson bra-kets. 



g(i,x) = (r---(-i,x)|^ar™*?5|rx'(+i,x)> 

= tr(r°"''^^75,S^(-i,x)7^r™*,S(i,x)r^"'=/i'^(-i,f,x)) (15) 
Here we measure the mass distribution for a static vector meson: 

python qcdutils_run . py \ 

-gauge : start = load : load = deino/demo . nersc . mdp \ 
-quark : kappa=0 .135: source_point=center: load=true \ 

-current_static : source_gamina = l : sink_gamina = l : curr ent_gamma = I : vtk=true \ 

> run, log 

Using the same diagram we can compute the spatial distribution of B* — )■ Bn by inserting 
the axial current {q^^q) in between a static B (q^^h) and and a static B* (q^y^h): 

python qcdutils_run . py \ 

-gauge : start=load : load=demo/ demo . nersc . mdp \ 
-quark : kappa=0 .135: source_point=center: load=true \ 

-current_static : source_gamma=l : sink_gamma=5 : curr ent_ gamma =5 : vtk=true \ 

> run, log nil,, —I ■■imia - ■ -I III. ■■■■Ill ■■Mil I— la ■■—■■nil iii ■iimiiiimiim ■■■■■■■■■■iiiii — ,„_ 1,11 

A sample image is shown in fig. 10. 
5.10 Localized instantons 

FermiQCD allows the creation of custom gauge configurations with localized topological 
charge. Here we consider the case of a pion propagator on a single gauge configuration in 
presence of one t'Hooft instanton (localized lump of topological charge). Here is the code: 
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Figure 10: Contour plot for a three points correlation function. 



python qcdutils_run . py \ 




-gauge : start=instantons : nt=20 : nx=20 : tO 


=0 : x0=4 . 5 : yO=10 : zO=10 \ 


-topcharge_vtk \ 




-quark :kappa=0. 120: source_point=center 


\ 


-pion : vtk=true > run . log 





This code places the center of the instanton at point (to, Xq, yo, zq) = (0, 4.5, 10, 10) and then 
computes a pion propagator with source on time slice but spatial coordinates {x,y,z) = 
(10,10,10) (center). 

Fig. 11 show the pion propagator in presence of the instanton as the instanton nears the 
center of the propagator. Each image has been generated using the above command by 
placing the instanton at different locations. The last image shows a superposition of the 
pion propagator with and without the instanton in order to emphasize the difference. The 
difference is small but visible. The propagator retracts as the instanton nears. One may say 
that the quark interacts with the instanton and acquires mass thus making the propagator 
decrease faster when going through the instanton. Fig. 12 shows the effect of the instanton 
on individual components of the quark propagator. 
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Figure 11: Pion propagator on a semi-cold configuration in presence of one localized in- 
stanton. The bottom right image shows the instanton (in blue) and an overlay of the pion 
propagator with and without the instanton. The difference shows that propagator shrinks as 
the instanton nears. 

6 Analysis with qcdutils.boot .py, qcdutils_plot .py, qcdutils.f it .py 

The console output of the qcdutils_run program consists of human readble text with com- 
ments and results of measurments performed on each gauge configuration. Here are some 
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Figure 12: Comparion of quark propagator components on a cold configuration (left) and in 
presense of a localized instanton (right). The instanton is located as in fig. 11 (bottom-right). 

examples of measurements logged in the output: 

plaquette = 0.654346 
C2 [0] = 14.5234 
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5 



4 



C3 [0] [0] = 1 .214321 
C3 [0] [0] = 1 . 123425 



6 



qcdutils_boot .py is a tool that can extract the vahies for these measurements, aggregate 
them and analyze them in various ways. For example it computes the average and bootstrap 
errors [16] of any function of the measurements. qcdutils_plot.py is a tool to visualize the 
results of the analysis. It uses the Python matplotlib package, one of the most powerful and 
versatile plotting libraries available, although qcdutils.plot .py uses only a small subset of 
the available functionality. 

Now, let us consider a typical Lattice QCD computation where one or more observables are 
measured on each Markov Chain Monte Carlo step (on each gauge configuration). We label 
the observables with Yj (gauge configuration, 2-point correlation function for a given value 
of t, etc.) 

We also refer to each measurements with yij where i labels the gauge configuration and j 
labels the observable (same index as Yj). yiQ could be, for example, the plaquette on the i th 
gauge configuration. 

The expectation value of each one observable is computed by averaging its measurements 
over the MCMC steps: 



Here N is the number of the measurements. The statistical error on the average for this 
simple case can be estimated using the following formula: 



Usually we are interested in the expectation value of non-trivial functions of the observables: 



Often the yij are not normal distributed and may depend on each other therefore standard 
error analysis does not apply. 




(16) 




(17) 



/ = (0| /(Fi, F2, ...Tm) |0) = /(Fi, F2, ...,1m) 



(18) 
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The proper technique for estimating the error on / is the bootstrap algorithm. It consists of 
the following steps: 

• We build K vectors of size N. The elements of these vectors are chosen at random, 
uniformly between {1,2, ...iV}. 

• For every k we compute: 

^E%. (19) 

i 

• Again for each k we compute: 

f^^f{Y^,Yl...,Y^) (20) 

• We then sort the resulting values for 

• We define the a percent confidence interval as [f^\ f^"\ where k' — [(1 — a)K/2\ and 
k" = \ l + a)K/2\. 

qcdutils_boot is a program that takes as input /(lo,^!,--) in the form of a mathemat- 
ical expression where the Yj are represented by their string pattern. It locates and ex- 
tracts the corresponding ijij values from the log files and stores them in a file called "qc- 
dutils_raw.csv" . It computes the autocorrelation for each of the and stores them in 
"qcdutils -autocorrelation. csv". It compute the moving averages for each of the Yj and stores 
them in "qcdutils_trails.csv". It generates the K bootstrap samples and saves them in 
"qcdutils_samplcs.csv". Finally it compute the mean and the 68% confidence level intervals 
[f^ 1 ] stores it "qcdutilsjresults.csv". 

Mind that these files arc created in the current working directory and they arc overwritten 
every time the qcdutils_boot is run. Move them somewhere else to preserve them. 

Moreover, if the input expression for / depends on wildcards, the program repeats the analysis 
for all matching expressions. 

qcdutils.boot performs this analysis without need to write any code. It only needs the input 
/ in the syntax explained below and the list of log-files to analyze for data. 
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6.1 A simple example 

Consider the output of one of the previous qcdutils_run: 

python qcdutils_run . py -gauge : load=* . mdp -plaquette > run . log 

In this case the observable is yo="plaquette". We can analyze it with 



This produces the following output: 

< plaquette > = min : 0.26, mean: 0.32, max: 0.35 
average trails saved in qcdutils_trails . csv 

bootstrap samples saved in qcdutils_samples . csv 
results saved in qcdutils_results . csv 



Notice that qcdutils_run takes three arguments: 

• A file name or file pattern (for example "run.log") 

• An expression (for example "plaquette"). 

• A condition (optional) 

Each of the argument must be enclosed in single quotes. 

The represents /(lo, ^i, ■■■) and the Yj are the names of observables in double quotes. 

In '"plaquette"' the outer single quote delimits the expression and the term plaquette 
between double quotes, determines the string we want to parse from the in file. 

qcdutils uses the observable name to find all the occurrences of 
plaquette = ... 



or 



in the input files and maps them into y^o where i labels the occurrence. In this case we have 
a single observable (plaquette) so we use to label it. 

The program opens the file or the files matching the file patterns and parses them for the 
values of the "plaquette" thus filling an internal table of y s. It gives the output as the result: 
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Here "mean" is the mean of the expression "plaquette" . min and max are the 65% confidence 
intervals computed using the bootstrap. 

Here is example of the content of the "qcdutils_results.csv" file for the average plaquette case: 

1 "plaquette" , " [min] " , " [mean] ", " [max] " 

2 "ylaquette " . . 26 . . 32 , . 38 

In general it contains one row for each matching expression. 

You can plot the content of the files generated by qcdutils_boot using qcdutils_plot: 



Here -r indicates that we want to plot the raw data, -a indicates we want a plot of auto- 
correlations, -t is for partial averages, and -b means we want a plot of bootstrap samples. 
qcdutils_plot loops over all the files reads the data in them and for each Yj it makes one 
plot with raw data (yij), one with autocorrelations, one with partial averages. Then for each 
/ it makes one plot with the bootstrap samples, and one plot with the final results found in 
"qcdutils_results.csv" . 

The plots are in PNG files which have a name prefix equal to the name of the data source 
file, followed by a serialization of the expression for Yj or /, depending on the case. 

For example in the case of the plaquette, the autocorrelations and the partial averages are 
in the files: 

1 qcdutils_autocorrelations_plaquette . png 

2 qcdutils_trails_plaquette . png 



and they are shown in fig. 13. 

Similarly, if you want to bootstrap /(Iq) 
run: 



= exp{YQ/3) where Yq is the plaquette you would 



1 ffi^KbJt^ □cdutiIs_bQOt ■ Dv^aaya,.^log. '.exn C ".vlaa.uette "/3l ' 

It produces output like this: 

Notice that again the observable Yq is identified for convenience by "plaquette". The double 
quotes are necessary to avoid naming conflicts between patterns and functions. 

Also notice that running qcdutils.boot twice does not guarantee generating the same exact 
results twice. That is because the bootstrap samples are random. 
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Figure 13: Example plots of autocorrelation (left) and partial averages (right). 
6.2 2-point and 3-point correlation functions 

In order to explain more complex cases we could generate 2- and 3- points correlation func- 
tions using something like: 



1 


python 


qcdutils_run . py \ 






2 




-gauge : start=cold : beta 


=4 


n=10 : steps=5 : therm=100 \ 


3 




-quark : kappa=0 . 11 : c_sw 


= 


4:save=false -pion 


4 




-4 quark : operator=5Ix5I 


> 


run . log 



For testing purposes can also run: 
1 python qcdutils_boot . py -t 



Where -t stands for test. This creates and analyzes a file called "test_samples.log" which 
contains random measurements for C2 and C3. Once this file is being created we can filter 
and study, for example, only the 2-point correlation function C2: 

1 python qcdutils_boot . py run . log '"C2[<t>]"' 

Notice that <t> means we wish to define a variable t to be used internally for the analysis and 
whose values are to be determined by pattern-matching the data. The t correspond to the 
j of the previous abstract discussion. "C2[<t>]" matches C2[0] with t = 0, C2[l] matches 
with t = 1, etc. 

The command above produces something like: 
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1 reading file test_samples.log 

2 C2 [00] occurs 100 times 

3 . . . 

4 C2 [15] occurs 100 times 

5 raw data saved in qcdutils_raw_data . csv 

6 autocorrelation for C2 [02] and d=l is -0.180453 

7 . . . 

8 autocorrelation for C2 [06] and d = l is -0.0436378 

9 autocorrelations saved in qcdutils_autocorrelations . csv 



10 


< 


C2 [00] 


> 




min : 


1 


988 , 


mean : 


1 


. 999 , 


max 


2 


. 008 


11 


< 


C2 [01] 


> 




min : 


1 


617 , 


mean : 


1 


.629 , 


max 


1 


. 64 


12 


< 


C2 [02] 


> 




min : 


1 


328 , 


mean : 


1 


.345 , 


max 


1 


. 359 


13 


< 


C2 [03] 


> 




min : 


1 


064 , 


mean : 


1 


.079 , 


max 


1 


.094 


14 


< 


C2 [04] 


> 




min : 





878 , 


mean : 





.894 , 


max 





. 908 


15 


< 


C2 [05] 


> 




min : 





722 , 


mean : 





.733 , 


max 





. 744 


16 


< 


C2 [06] 


> 




min : 





574 , 


mean : 





.584 , 


max 





.597 


17 


< 


C2 [07] 


> 




min : 





478 , 


mean : 





.49, 


max : 


0. 


5 


IS 


< 


C2 [08] 


> 




min : 





395 , 


mean : 





.407 , 


max 





.419 


19 


< 


C2 [09] 


> 




min : 





322 , 


mean : 





.331 , 


max 





.339 


20 


< 


C2 [10] 


> 




min : 





268 , 


mean : 





.277 , 


max 





.286 


21 


< 


C2 [11] 


> 




min : 





225 , 


mean : 





.231 , 


max 





.237 


22 


< 


C2 [12] 


> 




min : 





18 , mean : 


0. 


186 , 


max : 


0. 


192 


23 


< 


C2 [13] 


> 




min : 





138 , 


mean : 





. 144 , 


max 





. 151 


24 


< 


C2 [14] 


> 




min : 





107 , 


mean : 





. 112 , 


max 





. 118 


25 


< 


C2 [15] 


> 




min : 





0883 , 


mean 




0.0933, max: 


0.0 



26 average trails saved in qcdutils_trails . csv 

27 bootstrap samples saved in qcdutils_samples . csv 

28 ^m.^^2.t^-..^ved in acdutila-.xesull^s.».iasaa..^— 



which we can plot as usual with 



This produces about 60 plots. Some of them are shown in fig. 14. 

We can as easily compute the log of C2 (for every t): 
1 python qcdutils_boot . py test_samples.log ' log ( " C2 [<t >] " ) [_ 



or the log of the ratio between C2 at two consecutive time-slices: 

1 python qcdut ils_boot . py run2.log \ 

2 • log ("C2[<tl>] "/"C2[<t2>] ") • \ 

3 •t2==tl+l if tl<8 else t2==tl-l' 



In this case we used two implicit variables tl and t2 but we used the third argument of 
qcdutils_boot to set a condition to link the two. This produces the following output: 
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Figure 14: Example plots for the raw data (left), the distribution of raw data (center) and 
autocorrelations (right) for C2. 

1 reading file test_samples.log 

2 C2 [00] occurs 200 times 

4 C2 [15] occurs 200 times 

5 raw data saved in qcdutils_raw_data . csv 

6 autocorrelation for C2 [02] and d = l is -0.176706 

7 ... 

8 autocorrelation for C2 [06] and d = l is -0.0476376 

9 autocorrelations saved in qcdutils_autocorrelations . csv 

10 < log (C2 [00] /C2 [01] ) > = min : 0.196, mean: 0.204, max: 0.212 

11 < log (C2 [01] /C2 [02] ) > = min: 0.18, mean: 0.191, max: 0.202 

12 [ . . . ] 

13 < log (C2 [13] /C2 [14] ) > = min: 0.208, mean: 0.249, max: 0.291 

14 < log(C2 [14] /C2 [15] ) > = min: 0.147, mean: 0.189, max: 0.227 

15 average trails saved in qcdutils_trails . csv 

16 bootstrap samples saved in qcdutils_samples . csv 

17 results saved in qcdutils_results . csv 



In the same fashion we can compute a matrix element as the ratio between a 3-point corre- 
lation function (C3) and a 2-point correlation function (C2): 

1 python qcdutils_boot . py test_samples.log \ 

2 ' "C3[<t>] [<tl>]"/"C2[<t2>]"/"C2[<t3>]" ' \ 

3 ' t3==t and t2==t and tl==t ' > run. log 

4 python qcdutils_plot . py -a -t -b -r 



Some of the generated plots can be seen in fig. 19-16. 
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Figure 15: Example plots of moving averages (left) and distribution of bootstrap samples 
(right) for the ratio C3/C22. 
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Figure 16: Example plot showing results of the bootstrap analysis. 
6.3 Fitting data with qcdutiis_fit.py 

qcdutils_f it .py is a fitting and extrapolation utility. It can read and understand the output 
of qcdutils_boot .py. Internally it uses a "stabilized" multidimensional Newton method to 
minimize x^- It is stabilized by reverting to the steepest descent in case the Newton step 
fails to reduce the . The length of the steepest descent step is adjusted dynamically to 
guarantee that each step of the algorithm reduces the x^- The program accepts for input any 
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function and any number of the parameters. It also accepts, optionally, Bayesian priors for 
those parameters and they can be used to further stabilize the fit [17]. A more sophisticated 
approach is described in ref. [18]. 

In Euclidean space C2 can be modeled by an exponential aexp{—bt) and b is the mass of the 
lowest energy state which propagates between the source and the sink. Here is an example 
in which we fit C2 using a single exponential: 
python qcdutils_boot . py -t 

python qcdutils_boot . py test_samples.log '"C2[<t>]"' > run. log 
python qcdutils_f it . py ' a* exp ( -b *t) @a=2 , b =0 . 3 ' 

The input data is read from the output of qcdutils_boot. The expression in quotes is the 
fitting formula. You can name the fitting parameters as you wish (in this case a and b) but 
the other parameters (in this case t) must match the parameters defined in the argument 
of qcdutils_boot (<t>). The symbol separates the fitting function (left) from the initial 
estimates for the fitting parameters (on the right, separated by commas). Every parameter 
to be determined by the fit must have an initial value. 

The output looks something like this: 

a = 1.99864 

b = 0.200645 

Chi2= 12.8048378376 




qcdutils_f it .py also generates the plot of fig. 17 (left). 

If C2 is a meson propagator, b here represents the mass of the meson (of the lowest energy 
state with the same quantum numbers as the operator used to create the meson). 

Similary we can analyze and fit the log of C2: 

python qcdutils_boot . py test_samples.log ' log ( " C2 [<t >] " ) ' > run . log 
python qcdutils_f it . py ' a-b*t@a=l ,b=0 . 3 ' 

which produces something like: 

a = . 69169 

b = 0.200627 

chi2= 12.1641201448 

^iliTiiiiiiiliMmM^^ 

and the plot of fig. 17 (right) 

If our goal is obtaining b we can also cancel the a dependency in the analysis: 
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Figure 17: Example fits for a two points correlation function (left) and its log (right). 



python qcdutils_boot . py test_samples.log 


\ 


' log ("C2[<t>] "/"C2[<tl>] ") ' 'tl== 


t+1 ' > run . log 


python qcdut il s_f it . py 'b@b=0' 




and obtain: 


b = 0.201755 




chi2= 12.4502446913 




chi2/dof= 0.9577111301 





The generated plot is shown in fig. 18. 

Notice that the variable names a and b are arbitrary and you can choose any name. 
Similarly we can fit 3-point correlation functions: 

python qcdutils_boot . py test_samples.log ' " C3 [<tl >] [<t2>] " ' > run . log 
python qcdutils_f it . py ' a*exp (-b* (tl+t2) ) @a=3 , b=0 . 3 , _b=0 . 2 ' 

In this case we have stabilized the plot with a Bayesian prior, indicated by _b. A variable 
starting with underscore indicates the uncertainty associated with our a priori knowledge 
about the corresponding variable without underscore. In other words b=0.3,_b=0.2 is equiv- 
alent to b=0.3 ± 0.2. The result of this fit yields something like: 

a = 3.78387 

b = 0.195542 

chi2= 2070.73759118 

chi2/dof= 8.18473356199 
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Figure 18: Example plot of fit of log{C2{t) /C2{t - 1)). 
A call to qcdutils_plot .py generates the plot of fig. 19 (left) 

In order to extract a matrix element (for example a 4-quark operator) we fit the ratio between 
C3 and C2: 



python qcdutils_boot . py test_samples 


•log \ 


' "C3 [<t>] [<tl >] "/"C2[<t2>] "/" 


C2[<t3>] " • \ 


't3==t and t2==t and tl==t ' > 


run . log 


python qcdut il s_f it . py ' a@a=0 ' 




It produces output like: 


a = 1.00658 




chi2= 13.2075581022 




chi2/dof= 0.943397007297 





It produces the plot in fig. 19 (right). 

You can use qcdutils_f it .py to perform exatrapolations by using the -extrapolate command 
line option: 

python qcdutils_f it . py -extrapolate x = 100 ' ax + b@a=l , h =0 ' 

The extrapolated point will be added to the generated plot and represented by a square. 
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Figure 19: Example plot showing C3[t][tl] (left) and the fit of C3[t] [tl]/C2[t]c2[tl] (right). 
6.4 Dimensional analysis and error propagation 

In this section we did not discuss error propagation but we have developed a utility called 
Buckingham which is avalable from: 

http : / /code . google . com/p/buckingham/ 

It provides dimensional analysis, unit conversion, and aritmetic operation with error propga- 
tion. We plan to discuss it in a separate manual but we here provide one example of usage 
(from inside a Python shell): 

1 >>> from buckingham import * 

2 >>> a = Number(2.0, error=0.3, dims = "/ermi ") 

3 >>> b = NumberCl.O, error=0.2, dims= "second "2 ") 

4 >>> c = a/b 

5 >>> print c, c. units () 

(i (2.000 pm 0.500)/10"15 meter*second"-2 

7 >>> print c.convertC ' f ermi* second' -2 ') 

8 2 . 000 pm . 500 

9 >>> print C.convertC ' I i ghty ear * day ~ -2 ' ) 
10 (1 . 578 pm . 395) /10-21 

(here pm stands for ±) Buckingham supports 944 unit types (including eV) and their combi- 
nations. 
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A Filename conventions 



1 


Gauge configuration in NERSC format (3x3 or 3x2) 




2 


* . nersc 




3 


Gauge configuration in Fermiqcd format 




4 


* . mdp 




5 


Gauge configuration in MILC format 




6 


* . mile 




7 
8 


Generic LIME file 
* . lime 




9 


Gauge configuration in ILDG format 




10 


* . ildg 




11 


SciDAC quark propagator 




12 


* . s c idac 




13 


Quark propagator in FermiQCD format 




14 


* . prop . mdp 




15 


Time slice for gauge configuration in FermiQCD format 




16 


* . t [NMNN] . mdp 




17 


Time slice for propagator in FermiQCD format: 




18 


* . t [NNNN] . prop . mdp 




19 


Quark field for a given SPIN, COLOR source: 




20 


*. s [SPIN] . c [COLOR] . quark 




21 


Generic log file 




22 


* • log 




23 


VTK file containing real trace of plaquettes 




24 


* . plaquette . vtk 




25 


VTK file containing real part of Polyakov lines 




20 


* . polyakov . vtk 




27 


VTK file containing topological charge density 




28 


* . topcharge . vtk 




29 


VTK file containing topological charge density for a 


cooled config 


30 


* . topcharge . cool [STEP] . vtk 




31 


VTK file containt a the norm squared of a pion propagator 


32 


* . pion . vtk 




33 


HTML file generated by qcdutils_vtk , represents a VTK 


file . 


34 


* . vtk . html 




35 


Visit visualization script generated by qcdutils_vis . 


py 


36 


qcdutils_vis_ [UUID] . py 




37 


Visit image generates by the previous script 




38 


qcdutils_vis_ [UUID] _ [FRAME] . jpeg 




39 


Raw data extract from a log file by qcdutils_boot 




40 


qcdutils_raw_data . csv 




41 


Autocorrelations computed by qcdutils_boot 




42 


qcdutils_autocorrelations . csv 




43 


Partial averages computed by qcdutils_boot 
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44 qcdutils_trails . CSV 

45 Bootstrap samples generated by qcdutils_boot 

46 qcdutils_samples . CSV 

47 Means and bootstrap errors computed by qcdutils_boot 

48 qcdutils_results . CSV 



B Help Pages 

B.l qcdutils_get . py 



10 

11 

12 
13 
14 
15 
16 
17 
18 
19 
20 

21 
22 
23 
24 



$ qcdutils_get . py -h 






Usage : 








qcdutils_get 


■py 


[options] 


sources 


Examples : 








qcdutils_get 


■py 


--test 




qcdutils_get 


■ py 


-- convert 


ildg gauge . cold . 12x8x8x8 


qcdutils_get 


■py 


--convert 


mdp --float *.ildg 


qcdutils_get 


■py 


--convert 


split . mdp * . mdp 


Options : 








-h, --help 




show 


this help message and exit 


-q, --quiet 




no progress bars 


-d DESTINATION 




destination=DESTINATION 






destination folder 


-c CONVERT, -- 


convert=CDNVERT 






converts a field to format 






(ildg 


, split . prop . mdp , prop . ildg , prop . mdp , split . mdp , 






mdp ) 


-4, — float 




converts to float precision 


-8, --double 




converts to double precision 


-t , --tests 




runs 


some tests 


-n, --noprogre 


ssbar disable progress bar 



B.2 qcdutils_run.py 



1 $ qcdutils_run . py -h 

2 qcdutils_run . py is a tool to help you download and use fermiqcd from 

3 

4 http://code.google.eom/p/fermiqcd 

5 
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6 When you run : 

7 

8 python qcdutils_run . py [args] 

9 

10 It will : 

11 - create a folder called fermiqcd/ in the current working directory 

12 - connect to google code and download f ermiqcd . cpp + required libraries 

13 - if -mpi in [args] compile fermiqcd with mpiCC else with g++ 

14 - if -mpi in [args] run fermiqcd.exe with mpiCC else run it normally 

15 - pass the [args] to the compiled fermiqcd.exe 

16 

17 Some [args] are handled by qcdutils_run . py : 

18 -download force downloading of the libraries 

19 -compile force recompiling of code 

20 -source runs and compiles a different source file 

21 -mpi for use with mpi (mpiCC and mpirun but be installed) 

22 

23 Other [args] are handled by fermiqcd. cpp for example 

24 -cold make a cold gauge configuration 

25 -load load a gauge configuration 

26 -quark make a quark 

27 -pion make a pion 

28 (run it with no options for a longer list of options) 

29 

30 You can find the source code in fermiqcd/f ermiqcd . cpp 

31 

32 More examples : 



33 


qcdutils 


_run 


py 


-gauge 


start 


= cold 


: nt=16 : nx=4 




34 


qcdutils 


_run 


py 


-gauge 


start 


= hot : 


nt=16 : nx=4 




35 


qcdutils 


_run 


py 


-gauge 


load = 


cold . 


mdp 




36 


qcdutils 


_run 


py 


-gauge 


load = 


cold . 


mdp : steps=10: beta=5.7 


37 


qcdutils 


_run 


py 


-gauge 


load = 


* . mdp 


-plaquette 




38 


qcdutils 


_run 


py 


-gauge 


load = 


* . mdp 


-plaquette_vtk 




39 


qcdutils 


_run 


py 


-gauge 


load = 


* . mdp 


-polyakov_vtk 




40 


qcdutils 


_run 


py 


-gauge 


load = 


* . mdp 


- cool : steps =20 


-topcharge_vtk 


41 


qcdutils 


_run 


py 


-gauge 


load = 


* . mdp 


-quark : kappa=0 


12 : alg=minres_vtk 


42 


qcdutils 


_run 


py 


-gauge 


load = 


* . mdp 


-quark : kappa=0 


12 -pion 


43 


qcdutils 


_run 


py 


-gauge 


load = 


* . mdp 


-quark : kappa=0 


12 -pion_vtk 



44 

45 Options : 

46 -cool 

47 alg = ape 

48 alpha = 0.7 

49 steps = 20 
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-cool_vtk 




n = 20 




alpha = 0.7 




steps = 1 




cooling = 10 




-quark 




action = clover_fast (default) or clover_slow or clover_sse2 


alg = bicgstab (default) or minres or bicgstab_vtk or 


minres_vtk 


abs_precision = le-12 




rel_precision = le-12 




source_t = 




source_x = 




source_y = 




source_z = 




source_point = zero (default) or center 




load = false (default) or true 




save = true (default) or false 




matrices = FERMILAB (default) or MILC or 




UKQCD or Minkowsy -Dirac or Minkowsy - Chiral 




kappa = 0.12 




kappa_t = quark [ "fcaypa "] 




kappa_s = quark [ "fcappa "] 




r_t = 1.0 




r_s = 1.0 




c_sw = 0.0 




c_E = 0.0 




c_B = 0.0 




-meson 




source = 5 




sink = 5 




current = I 




-4quark 




source = 5 (default) or I or or 1 or 2 or 3 or 05 or 




15 or 25 or 35 or 01 or 02 or 03 or 12 or 13 


or 23 


operator = 51x51 (default) or 01x01 or 11x11 or 21x21 


or 31x31 or 


051x051 or 151x151 or 251x251 or 351x351 or 


011x011 or 


021x021 or 031x031 or 121x121 or 131x131 or 


231x231 or 


5Tx5T or OTxOT or ITxlT or 2Tx2T or 3Tx3T or 05Tx05T or 


15Txl5T or 25Tx25T or 35Tx35T or OlTxOlT or 


02Tx02T or 


03Tx03T or 12Txl2T or 13Txl3T or 23Tx23T 




-gauge 




nt = 16 




nx = 4 




ny = nx 
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96 




start = load (default) or cold or hot or instantons 


97 




load = demo.mdp 


98 




n = 


99 




steps = 1 


100 




therm = 10 


101 




beta = 


102 




zeta = 1.0 


103 




u_t = 1.0 


104 




u_s = 1.0 


105 




prefix = 


106 




action = wilson (default) or wilson_improved or wilson_sse2 


107 




save = true 


108 




to = 


109 




xO = 


110 




yo = 


111 




zO = 


112 




rO = 1.0 


113 




tl = 1 


114 




xl = 1 


116 




yl = 1 


116 




zl = 1 


117 




rl = 0.0 


118 




-baryon 


119 




-pion 


120 




-pion_vtk 


121 




-meson_vtk 


122 




-current _static 


123 




-current_static_vtk 


124 




-plaquette 


125 




— plaquette vtk 


126 




-polyakov_vtk 






-topcharge_vtk 




B.3 


qcdutils_vis .py 


1 


$ qcdutils_vis . py -h 


2 


Usag 


e : 


3 


This 


is a utility script to manipulate vtk files containing scalar files. 


4 


Files can be split , interpolated , and converted to jpeg images . 


5 


The 


conversion to jpeg is done by dynamically generating a visit script 



6 that reads the files, and computes optimal contour plots. 

7 

8 Examples : 

9 

10 1) make a dummy vtk file 
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11 

12 qcdutils_vis . py -m 10 f older / test . vtk 

13 

14 2) reads fields from multiple vtk files 

16 

16 qcdutils_vis . py -r field folder/*. vtk 

17 

18 3) extract fields as multiple files 

19 

20 qcdutils_vis . py -s field folder/*. vtk 

21 

22 (fields in files will be renamed as "slice") 

23 4) interpolate vtk files 

24 



26 qcdutils_vis . py -i 9 folder/*. vtk 

26 

27 tricubic Resample / Interpolate individual vtk files 

28 

29 visit -V 10x10x10 folder/*. vtk 

30 

31 6) render a vtk file as a jpeg image 

32 

33 qcdutils_vis . py -p ' AnnotationAttributes [axes3D . bboxFlag =0] ; 

34 Res amp leAttrihutes [s amp lesX=160 ; s amp lesY=160; s amp lesZ = 160] ; 

36 ContourAttributes [SetMultiColor (9 , $orange)] ' ' folder/* . vtk ' 

36 

37 or simply 

38 

39 qcdutils_vis . py -p default ' folder/* . vtk ' 

40 

41 Options : 

42 -h, --help show this help message and exit 

43 -r READ, --read=READ name of the field to read from the vtk file 

44 -s SPLIT, — split = SPLIT 

46 name of the field to split from the vtk file 

46 -i INTERPOLATE, — interpolate = INTERPOLATE 

47 name of the vtk files to add/interpolate 

48 -c CUBIC, --cubic-interpolate=CUBIC 

49 new size for the lattice 10x10x10 

50 -m MAKE, --make=MAKE make a dummy vtk file with size"3 whete size if 

arg of 

51 make 

52 -p PIPELINE, --pipeline=PIPELINE 

53 visualizaiton pipeline instructions 
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B.4 qcdutils_vtk . py 



1 


$ qcdutils_vtk . py 


-h 






2 


Usage 


: qcdutils_vtk . py filename. vtk 






3 
4 


Options : 








5 


-h, 


--help 


show this help message and 


exit 


B 


-u 


UPPER, --upper -threshold=UPPER 






7 






treshold for isosurface 




8 


-1 


LOWER, --lower-threshold=LOWER 






9 






treshold for isosurface 




10 


-R 


UPPER_RED , - 


-upper -red=UPPER_RED 






11 






color component for 


upper 


isosurface 


12 


-G 


UPPER_GREEN , 


--upper -green=UPPER_GREEN 






13 






color component for 


upper 


isosurface 


14 


-B 


UPPER_BLUE , 


— upper -blue =UPPER_BLUE 






15 






color component for 


upper 


isosurface 


16 


-r 


LOWER_RED , - 


-lower -red=LOWER_RED 






17 






color component for 


lower 


isosurface 


18 


-g 


LOWER_GREEN , 


--lower -green=LOWER_GREEN 






19 






color component for 


lower 


isosurface 


20 


-b 


LOWER_BLUE , 


--lower- blue =LOWER_BLUE 






21 






color component for 


lower 


isosurface 



B.5 qcdutils_boot . py 



1 $ qcdutils_boot . py -h 

2 Usage: qcdutils_boot.py *.log ' x [<a>]/y [<h>] ' ' abs (a-b) ==1 ' 

3 scans all files *.log for expressions of the form 

4 X [<a>] =< value > and y [<b>] =<value > 

5 and computes the average and bootstrap errors of x [<a>] /y [<b>] 

6 where <a> and <b> satisfy the condition abs(a-b)==l. 

7 

8 This is program to scan the log files of a Markov Chain Monte Carlo 

Algorithm , 

9 parse for expressions and compute the average and bootstrap errors of any 

10 function of those expressions.lt also compute the convergence trails of 

the 

11 averages . 

12 

13 Options : 

14 --version show program's version number and exit 

15 -h, --help show this help message and exit 

16 . -b MIN. r.=.mjLja£Mmm^dP^.=MTN 
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17 






the first occurrence of expression to he 
consi dered 


18 


-e 


MAX, --maxmium_ 


index=MAX 


19 






the last occurrence +1 of expression to be 
cons i dered 


20 


-n 


NSAMPLES , — numh er_ of_ s amp I es=NSAMPLES 


21 






number of required bootstrap samples 


22 


-P 


PERCENT, — percentage=PERCENT 


23 






percentage in the lower and upper tails 


24 




--test 


make a test ! 


26 


-r. 


--raw 


Load raw data instead of parsing input 


26 


-a. 


--advanced 
variab I e 


In advanced mode use regular expressions for 


27 






patterns 


28 


-i 


IMPORT_MODULE , 


— import_module = IMPORT_MODULE 


29 






import a python module for expression evaluation 


30 
31 


-0 


OUTPUT_PREFIX , 


— output_prefix = OUTPUT_PREFIX 



B.6 qcdutils_plot . py 



1 $ qcdutils_plot . py -h 

2 Usage: python qcdutils_plot . py 

3 

4 plot the output of qcdutils.py 

5 

6 Options : 



7 --version show program's version number and exit 

8 -h , --help show this help message and exit 

9 -i INPUT_PREFIX , -- input _prefix=INPUT_PREFIX 

10 the prefix used to build input filenames 

11 -r , --raw make raw data plots 

12 -a, --autocorrelations 

13 make autocorrelation plots 

14 -t , --trails make trails plots 

15 -b , --bootstrap -samples 

16 make bootstrap samples plots 

17 -V PLOT_VARIABLES , — plot_variables=PLOT_VARIABLES 

18 plotting variables 

19 -R RANGE, --range=RANGE 

20 range as in 0:1000 ^^^^^^^^^^^^^^^^^ 



B.7 qcdutils_f it .py 
1 $ qcdutils_f it . py -h 



78 



2 Usage: qcdut i 1 s_f it . py [OPTIONS] ' expression@values ' 

3 Example: qcdutils -f it . py ' a*x + b@a=3 , b=0 ' 

4 default filename is qcdutils_results . csv 

5 . . . . , 'x ' , ' min ' , ' mean ' , 'max ' 

6 23, 10, 11, 12 

7 ..... etc etc etc 

8 

9 Options : 

10 

11 

12 
13 
14 
15 
16 
17 
18 
19 
20 
21 
22 
23 
24 
25 



--version show program's version numb ev and exzi 

-h , --help show this help message and exit 

-i INPUT, --input=INPUT 

input file (default qcdutils_results . csv) 
-c CONDITION, --condition=CONDITION 

sets a filter on the points to be fitted 
-p PLOT, --plot=PLOT plots the hessian (not implemented yet) 
-t , --test test a fit 

-e EXTRAPOLATIONS , --extrapolate=EXTRAPOLATIONS 

extrpolation point 
-a AP , --absolute_precision=AP 

absolute precision 
-r RP , --relative_precision=RP 

relative precision 
-n NS , --number_steps=NS 
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